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ABSTRACT 


The  energy  crisis  has  brought  about  public  concern  of 
the  limitation  of  conventional  resources.  Suggestions  have 
been  made  to  conserve  energy  and  to  seek  alternative 
resources.  Solar  energy  has  a  promising  future  to  serve  as 
alternative  resource  especially  in  the  area  of  space 
heat ing . 

The  controls  problem  of  a  solar  heating  system  is  shown 
to  have  a  significant  impact  on  energy  use.  A  new  approach 
to  this  problem  is  proposed  and  studied  in  this  thesis. 
Thermal  as  well  as  economic  performance  of  the  system  under 
this  new  suboptimal  controller  was  also  investigated. 
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1.  INTRODUCTION 


1.1  The  Importance  of  Energy. 

From  the  earliest  time, Man  has  tried  to  meet  his 
increasing  demands  for  goods  and  services  through  animal  and 
especially  mechanical  slaves;  tools,  machines  and  engines. 
Mechanical  slaves,  in  turn,  require  fuel  to  operate,  and 
fuel  appears  in  various  forms  such  as  coal,  oil,  uranium, 
etc  . 

Fuel  or  Energy  which  was  abundant  and  relatively 
inexpensive  is  fundamental  to  present  life  style;  our  food 
supply,  our  tremendous  industry,  our  conveniences  of  life, 
etc.,  ail  depend  upon  a  ready  supply  of  energy. 

1.2  The  Energy  crisis. 

The  heavy  dependence  of  the  industrialized  countries  on 
oil  has  made  them  vulnerable  to  any  happening  in  the  Middle 
East.  The  oil  embargoes  and  price  increases  of  the  oil 
cartel  OPEC  have  caused  chaos  and  panic  among  the  developed 
as  well  as  developing  countries.  That  has  ,in  turn, 
triggered  public  concern  about  the  limitation  of  the 
conventional  resources. 

To  solve  the  problem,  it  has  been  suggested  not  only 
that  more  domestic  explorations  into  traditional  and 
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alternative  sources  of  energy  be  promoted  but  that  the 
public  also  be  taught  the  importance  of  energy  conservation. 


Oilsand  projects  in  Alberta  have  appeared  quite 
promising  in  terms  of  making  Canada  independent  of  imported 
oil  for  at  least  into  the  next  century.  However,  this  would 
be  only  a  short  term  solution  to  the  "persistent"  energy 
problem,  for  oil  is  a  limited  source  of  energy. 

Is  coal  able  to  replace  oil  in  providing  us  with  an 
unaepletable  supply  of  energy?  The  answer  is  quite  clear. 
Although  coal  is  considered  one  of  the  most  abundant 
resources  in  North  America,  it  is  limited  in  quantity.  For 
example:  The  U.S.  has  an  estimated  3.2  trillion  tons  of  coal 
[ 1 ] ,  of  which  about  217  billion  are  economically  mineable. 
The  U.S.  will  run  out  of  an  economic  coal  supply  in  180 
years,  if  it  is  mined  at  the  rate  of  some  1200  million  tons 
per  year  as  proposed  by  the  Federal  Economic  Agency  (FEA). 

What  about  the  conservation  of  energy?  It's  not  a 
promising  solution  either,  for  it's  hard  to  convince  people 
to  give  up  their  life  styles  overnight.  Therefore,  new 
sources  of  energy  such  as  nuclear  power  and  solar  energy 
need  to  be  developed  to  replace  conventional  ones. 

Nuclear  power  can  certainly  satisfy  our  energy  demand 
but  there  are  still  doubts  about  its  safety.  The  Three-mile 
Island  accident  is  one  of  the  reasons  why  people  feel 
concerned.  Moreover,  it  also  depends  on  a  finite  resource, 


uranium. 
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On  the  contrary,  solar  energy  is  a  renewable  source  and 
poses  no  potential  danger  whatsoever  to  its  users.  Hence,  it 
seems  to  be  the  answer  to  our  persistent  problem. 


1.3  Solar  Energy,  a  new  source  of  energy. 

Looking  up  at  the  sky,  we  see  a  resource  that  can 
supply  us  with  all  the  energy  we  need:  the  inexhaustible, 
pol lut i on- f ree  energy  coming  from  the  sun.  How  extensive  is 
this  resource?  Lake  Erie  receives  enough  sunlight  to  provide 
the  total  U.S.  energy  supply. 

Before  a  new  source  of  energy  is  introduced  into  our 
traditional  way  of  life,  three  conditions  need  to 
considered;  namely,  it  should  be  technically  sound  and 
feasible,  economically  competitive  with  the  conventional 
sources  and  socially  accepted. 

Although  solar  energy  is  not  yet  the  answer  to  our 
energy  problem,  it  has  been  proven  technically  feasible  end 
economically  competitive  for  a  few  applications.  There  has 
been  increasing  interest  in  the  use  of  solar  energy  for 
heating  and  cooling,  since  the  amount  of  energy  saved  would 
be  considerable.  Space  heating  and  water  heating  account  for 
22%  of  U.S.  energy  consumption.  For  water  heating  alone,  the 
amount  saved  could  light  every  light  in  the  U.S.  twicefl]. 

Solar  energy  has  also  been  applied,  although  on  a 
limited  basis,  to  agriculture  and  industry[1].  It  has  been 
projected  that  about  5%  of  energy  demand  for  U.S. 
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agriculture  will  be  supplied  by  solar  energy  by  1985.  By 
2000,  the  estimated  percentage  will  be  25%,  and  by  2020  it 
will  become  50%. 

In  industrial  applications  such  as  industrial  hot  water 
and  drying/dehydration  systems,  it  is  estimated  that  about 
0.2%  of  the  energy  needed  for  such  applications  may  be 
supplied  by  solar  energy  by  1985,  or  about  20%  by  2020. 
Research  is  also  being  actively  pursued  into  electrical 
applications.  According  to  [8],  by  1983,  Bleyle  of  America 
Inc.,  will  receive  more  than  half  of  its  electricity  demands 
from  solar  energy,  and  this  project  is  expected  to  become 
cost-effective  by  some  time  between  1990  and  1995. 

Regarding  environmental  concerns,  solar  energy  has  not 
caused  any  such  serious  problems  as  nuclear  power,  thus 
making  it  socially  acceptable;  although  the  occupation  of 
large  areas  by  collectors  could  be  considered  detrimental  by 
some  owners. 


1.4  The  heating  problem. 

Space  heating  has  been  recognized  as  the  most 
favourable  area  of  solar  energy  applications.  In  a  study 
reported  by  Arnold  D.  Cohen  of  General  Electric[1],  16 
million  buildings  (21%  of  all  U.S.A.  buildings)  will  be 
equipped  with  solar  heating  systems  by  the  year  2000. 

It  is  generally  agreed  that  the  way  a  building  is 
operated  has  a  significant  impact  on  energy  use  and  that  a 
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poorly  operated  building  can  defeat  even  the  best 
energy-conserving  design.  It  is  also  accepted  that  the  use 
of  controls  will  lead  to  improved  system  efficiencies  and 
that  the  uses  of  these  controls  require  thorough 
i  nves t i ga t i on . 

The  approach  generally  used  to  solve  this  problem  is  to 
model  the  system  and  then  to  find  an  optimal  control 
strategy  for  the  model. 

The  active  solar  heating  systems  all  use  the  same  kinds 
of  components  although  they  might  have  many  different 
configurations.  A  solar  collector,  storage  tanks,  pumps,  a 
load  heat  exchanger,  and  an  auxiliary  heater  are  the 
components  that  seem  indispensable  in  an  active  solar 
heating  system. 

Each  component  has  a  different  function[2]. 

The  solar  collector  is  used  to  collect  energy  from  the 
sun  which  is  then  delivered  to  the  load  through  the  load 
heat  exchanger  or  stored  in  the  storage  tanks.  Supplemental 
energy  is  provided  by  the  auxiliary  heater  when  solar  energy 
is  depleted  or  insufficient.  Pumps  operate  to  carry  energy 
from  one  part  of  the  system  to  the  other.  The  successful 
operation  of  each  component  contributes  to  the  successful 
operation  of  the  whole  system.  Depending  on  the  desired 
level  of  complexity,  there  are  many  ways  to  model  those 
components.  According  to  D.M.  Auslander  et  al[2j: 

”...,the  load  is  often  modeled  as  a  single  'lump', 
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i.e.,  one  temperature  characterizes  the  instantaneous 
state  of  the  load.  The  storage  tank  is  also  modeled  as 
a  single  lump.  Varying  levels  of  complexity  can  be  used 
to  model  the  heat  exchangers,  we  have  chosen  to  use 
single  lump  models  for  each  of  the  heat  exchangers.  All 
flow  components(  pumps,  piping,  valves)  have  been 
modeled  as  static  elements.  The  collector  is  modeled  as 
a  heat  exchanger  in  which  the  heat  input  from 
insolation  is  a  decreasing  function  of  the  collector 
temperature.  Any  desired  function  can  be  used." 


The  controls  problem  is  considered  after  every  component  of 
the  system  has  been  modeled. 

The  controls  problem  is  concerned  with  how  to  control 
the  interaction  between  these  components  to  lead  to  a 
successful  design.  It  is  generally  agreed  that  the  controls 
problem  consists  of  a  regulation  problem  and  a  minimization 
problem [ 2 ] . 

Auslander  has  argued  that  the  primary  obstacle  to 
developing  a  method  of  comparing  competing  controller 
designs  is  due  to  the  lack  of  genera]  agreement  on  a  uniform 
standard  or  specifications. 

After  rejecting  as  unreasonable  the  quadratic 
performance  index  which  combines  temperature  deviation  from 
the  set  point  and  other  factors  such  as  auxiliary  energy 
usage,  Auslander  went  on: 


"Rather  than  a  quadratic  performance  index,  then,  it 
seems  reasonable  to  base  controller  comparison  on  some 
temperature  band  or  minimum  temperature...  Using  the 
minimum  temperature  as  standard,  two  competing 
controllers  could  be  compared  by  adjusting  them  so  that 
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they  both  have  the  same  minimum  temperature  over  the 
period  of  interest.  Auxiliary  energy  usage,  or  other 
performance  factors,  could  then  be  compared 
directly ..." 


A  good  controller  is  one  that  not  only  regulates  the  system 
with  respect  to  disturbances,  solar  heat  input  and  ambient 
temperature,  but  also  conforms  to  the  minimization 
standards.  According  to  Winn[3],  there  are  three  levels  of 
minimization. 

His  first  level  is  minimization  of  the  total  heat  input 
into  the  load  (not  just  the  auxiliary  heat).  The  second 
level  is  minimization  of  the  auxiliary  heat  used  while 
meeting  the  regulation  spec i f icat ions .  The  last  level  is 
minimization  of  parasitic  energy. 

The  second  level  of  minimization  is  the  one  that  will 
be  used  as  as  minimization  standard  in  this  thesis.  Reasons 
for  not  using  the  first  and  third  levels  have  been  explained 
in  detail  by  Au slander  et  al  in  [2]. 


1.5  Summary  of  past  results. 

In  [3]  E.  Winn  et  al.  have  divided  the  optimal  problem 
of  a  solar  energy  system  into  three  kinds: 


"1.  Optimal  controllers  of  the  First  Kind  are 
represented  by  those  controllers  that  optimally  supply 
heat  to  a  building  in  a  manner  such  that  a  measure  of 
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the  energy  supplied  and  the  occupant  discomfort  is 
minimized . 

2.  Optimal  controllers  of  the  Second  Kind  represent 
controllers  that  maximize  the  difference  between  the 
useful  energy  collected  and  the  pumping  costs 
associated  with  collecting  the  solar  energy. 

3. Optimal  controllers  of  the  Third  Kind  represent 
controllers  that  combine  the  collection  and 
distribution  function." 


Then,  Winn  proceeded  to  find  the  controllers  of  the  above 
three  kinds.  It  is  the  Second  Kind  controller  that  has 
practical  significance  for  the  collector-control  problem. 
According  to  Winn,  the  solution  to  the  Second  Kind  problem 
is  a  proportional  controller.  In  establishing  the  objective 
function  of  the  Second  Kind  optimal  controller,  Winn 
included  parasitic  losses  of  the  pumps;  however,  there  is 
some  doubt  that  pumping  power  should  be  considered  as  a 
loss,  since  with  proper  placement  of  components  this  power 
will  also  be  a  heat  input  to  the  system,  thus  offsetting  a 
corresponding  amount  of  auxiliary  heat  input[4]. 

There  have  been  debates  over  whether  or  not  the 
collector  pump  controller  should  be  on-off  or  proportional. 
Before  discussing  the  collector  pump  controller,  it  is 
necessary  to  describe  the  equation  of  the  collector.  A  flat 
plate  solar  collector  is  a  device  that  captures  heat  from 
the  sun,  the  useful  heat  delivered  by  a  solar  collector  is 
equal  to  the  energy  absorbed  by  the  collector  plate  less  the 
heat  lost  to  the  surroundings.  The  heat  rate  input  to  the 
system  from  the  collector  can  be  calculated  as  follows[9]: 
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Qc=Ac(Ht  ar  -UL(Tp-Ta ) )  (1.5-1) 

where,  Qc=useful  heat  delivered  by  the  collector 
Ac=  effective  collector  area  (m2) 

Ht=  incident  solar  radiation  (KJ/hr*m2) 

(received  on  the  tilted  collector  surface) 
t-  collector  cover  transmittance 
a  =  plate  absorptance 

U  =  collector  overall  loss  coef f ic i ent (KJ/hr°Cm2 ) 

XJ 

Tp =  average  absorber  plate  temperature 
Ta =  ambient  temperature 

The  term  A c H „ (  ar  )  is  the  amount  of  energy  absorbed  at 
the  absorber  surface.  This  quantity  depends  on  the  effective 
collector  area,  the  solar  radiation,  the  transmissivity  of 
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major  factors  that  change  the  transmissivity  coefficient.  It 

is  the  optical  property  of  the  collector  surface  that 

determines  the  absorptivity  coefficient.  Black  surfaces  have 

high  absorptivity  for  the  visible  range  of  the  solar 

spectrum.  The  absorptivity  coefficients  of  carbon  black, 

metal  oxides,  and  black  paints  often  have  values  above  0.95. 

The  heat  loss  from  the  collector  is  represented  by  the 

term  A  U  (T  -  T  ).The  overall  heat  loss  coefficient,  U  ,  is 
c  L  P  ^  t 

in  the  order  of  6  to  11  W/irrC  for  one  glass  glazing  and 
about  4  W/m2°C  for  two  panes. 
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In  equation  (1.5-1),  Tp  represents  the  average  plate 
temperature,  which  is  difficult  to  measure.  So,  instead  of 
using  Tp  ,  we  use  T  ,  which  is  the  temperature  of  the  fluid 
entering  the  collector,  to  calculate  the  heat  delivered  by 
the  collector. 

Qc=F  Ac(HTar-UL(T.-Ta ) )  (1.5-2) 

where,  F  is  a  correction  factor,  or  heat  recovery  factor. 
Its  value  is  between  0.  and  1.  such  that  Qc  calculated  by 
equation  (1.5-2)  is  equal  to  that  evaluated  by  equation 
(1.5-1). 


If  we  ignore  the  heat  loss  in  the  pipes  between  the 
collector  and  the  storage  tank,  i.e.  the  storage  tank  Ts 
equals  the  inlet  collector  temperature  Ti  ,  and  if  we 
incorporate  the  collector  pump  controller,  u3,  into  equation 


(1.5-2),  we  o brain: 

Qc=AcF  (u,;(HTar-UL(Ts-Ta)} 
-If  on-off  controller  is  used 


( 1 . 5- 3  ) 


F ( u 5 )  =  0 .  if  u  3  =  0  . 

F ( u  3 ) =F  if  u  3  =  u  3 max 

-If  proportional  controller  is  used,  value  of  F(u3)  is 
between  0.  and  F 

To  resolve  the  controversy  over  on-off  and  proportional 
controllers,  a  comparative  study  on  on-off  and  proportional 
controlled  system  has  been  reported  by  Lewis  and  Carr[5]. 

The  following  is  results  of  solar  collection  in  one  day 


of  the  two  controllers. 


1 1 


ou/qt 

qu/qt 

— T 

Proport i onal 

On-Off 

325  W/m2 

0.934 

0.955 

490. 

0.945 

0 .955 

620. 

0.947 

0 . 954 

70  1  . 

0.953 

0.954 

73  1  . 

0.955 

0.955 

where 


Thus  , 


Qu=  useful  energy  collected 

Qt=  total  energy  incident  on  absorber  plate. 

Ht=  solar  insolation  incident  on  absorber  plate, 
according  to  Lewis  et  al.: 


"The  slightly  lower  energy  gains  made  by  the 
proportionally  controlled  system  are  due  to  the  low 
fluid  flow  rates  through  the  system  at  low  temperature 
di  f  f erent  ials . . . 

While  this  low  flow  rate  proved  to  be  advantageous 
at  earlier  hours  by  enabling  the  system  to  gain  energy 
without  excessive  cycling,  it  is  now  of  some 
disadvantage  because  the  low  flow  rate  holds  plate 
temperature,  and,  therefore  heat  losses  higher  than 
those  encountered  in  the  on/off  system." 


McDonald  et  a  1 [ 6  ]  used  an  on-off  controller  on  the  collector 
side  and  used  an  adaptive  optimal  algorithm  to  control  other 
variables . 


The  approach  was  basically  to  identify 
model  of  the  system,  then  to  employ  optimal 


a  linearized 
control  theor v 


to  determine  gains  of  the  optimal  controller  which  minimize 


a  cost  function.  This  process  repeats  itself  in  the  next 
intervals  of  time.  As  McDonald  put  it: 
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"The  actual  building  and  HVAC  system  is  a  non-linear 
system  with  operating  points  which  can  vary  over  a  wide 
range.  The  linearized  model  is  valid  only  about  a 
region  of  the  operating  point;  thus,  the  identification 
of  a  linearized  mode  1  must  be  an  on-going  process  with 
optimal  controller  being  modified  or  adapted  for  each 
new  linearized  model  or  operating  point  of  the  system." 


The  cost  was  chosen  as  the  integral  quadratic  cost 
functional  of  plant  state  variables  and  control  variables. 
The  weighting  matrices  in  the  cost  functional  assign 
relative  importance  to  state  and  control  variables.  Large 
weights  were  given  to  room  temperature  and  auxiliary  energy 
variable  in  order  to  maintain  comfort  conditions  and  to 
minimize  auxiliary  energy  variables. 

A  conventional  approach  was  also  presented  to  compare 
its  result  with  that  of  the  Adaptive  Optimal  controller.  In 
Fia.1.1,  T,  represents  room  temperature  and  T  is  the  ambient 

b  3. 

temperature.  If  the  point  representing  room  and  ambient 
temperature  in  (T,  ,T  )  plane  is  in  the  hatched  region,  only 

fa 

solar  energy  will  be  used,  if  not  auxiliary  heat  will  also 
be  used.  Simulation  results  for  the  heating  season  are 
summarized  in  table  1.1. 

Room  temperature  was  shown  to  maintain  slightly  more  closely 
to  the  desired  by  the  conventional  controller  than  the 
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adaptive  one.  However,  the  adaptive  controller  used  28%  less 
auxiliary  energy  than  the  conventional  due  to  the  fact  that 
storage  tank  temperatures  were  kept  at  lower  values  ,  thus 
making  the  collector  operate  more  efficiently  and  collect 
more  solar  energy.  This  approach  has  once  more  solidified 
the  role  of  modern  control  theory  in  solar  energy  heating 
system  ;however,  it  requires  rather  sophisticated 
computations  at  each  update  and  is  inherently  limited  in  its 
optimization  look-ahead  time  because  of  the  limited  duration 
of  accuracy  for  the  linearized  model[4].  In  this  thesis  a 
more  direct  approach,  based  on  an  accurate  nonlinear  model, 
is  taken. 


1.6  Scope  of  this  thesis. 

A  sub-optimal  adaptive  controller  for  a  solar-assisted 
heat-pump  system  is  proposed.  Its  simulation  results  are 
shown  to  give  significant  improvement  over  those  of  WATSUN 
program  [7] (more  will  be  said  about  the  WATSUN  approach  in 
chapter  three).  In  view  of  the  results  cited  above,  the 
collector  pump  controller  is  taken  to  be  on-off.  Also, 
following  the  suggestion  of  Auslander  et  a 1 [ 2  ]  ,  the  building 
loop  controller  is  assumed  to  be  thermostatic.  The  major 
control  problem  remaining  is  then  just  the  control  of  the 
heat  pump  that  is  used  to  upgrade  the  collected  solar 


energy . 
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In  chapter  two,  the  system  models  are  developed  before 
the  optimization  problem  is  addressed.  Also  presented  in 
this  chapter  is  the  derivation  of  the  adaptive  preview 
suboptimal  controller.  The  stability  of  this  suboptimal 
controller  is  investigated  in  chapter  three.  The  simulation 
results  of  the  suboptimal  controller  and  the  WATSUN  approach 
are  compared  and  analysed  in  later  sections  of  chapter 
three.  In  the  last  chapter,  chapter  four,  the  thermal  as 
well  as  economic  performance  of  the  system  are  examined  when 
major  parameters  of  the  system  are  varied.  The  viability  of 
the  system  is  also  addressed  here. 

Although  the  approach  can  be  practically  adopted,  it  is 
advisable  to  make  a  much  more  thorough  investigation  before 
success  can  be  assured.  The  controller  has  been  designed 
with  the  assumption  that  the  weather  is  a  deterministic 
process,  or  at  least  the  average  ambient  temperature  of  the 
next  day  is  certainly  known.  However,  the  author  believes 
that  results  will  not  be  very  much  different  even  if  the 
next  day  average  ambient  temperature  cannot  be  predicted 
precisely.  In  other  words,  it  is  possible  to  have  similar 
results  even  if  the  process  is  stochastic  in  nature. 

This  controller  is  believed  to  be  suitable  for  eventual 
application  even  though  the  work  is  just  exploratory. 


15 


TABLE  1 .  1 
Heating  Energy 


(  10b 

/BTU) 

Aver,  room 

Aux  . 

solar 

I n ternal 

Total 

temp . 

%sav i ng 

CONV 

3.74 

1  .05 

6.12 

10.90 

70.01 

ADAP 

2.66 

2.06 

6.12 

10.84 

69.35 

28 . 8% 
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Fig. 1.1:  Conventional  Controller  schedule  for  solar  energy 


usage 
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2.  OPTIMIZATION 


2  .  1  I  n t roduct i on . 

Attention  has  recently  been  given  to  problems  of 
optimal  control  of  heat  storage  and  heat  transfer  processes 
in  HVAC  systems.  Control  objectives  are  to  minimize  the 
seasonal  consumption  of  purchased  energy  where  solar 
collection  and  storage  are  employed,  and  /or  to  minimize  the 
peak-period  consumption  of  electrical  energy  where  off-peak 
charging  are  utilized.  Among  the  typical  control  inputs  for, 
say,  a  solar  assisted  heat  pump  system  are  mass  flow  rate 
through  solar  collectors,  heat  pump  input  power,  flow  races 
to  building  convectors,  and  building  air-handling  variables. 
State  variables  are  typically  the  building  ana  storage 
temperatures . 

Potential  savings  with  optimal  control  seem  to  be 
considerable;  simulations  have  shown  that  improvements  of 
50%  or  mere  are  possible  under  certain  conditions  with 
optimal  control  of  solar  collector  and  building  variables, 
as  compared  with  conventional  cont rol [ 1 , 2 ] .  Better  design 
concepts  are  needed,  however,  for  the  realization  of 
practical,  near-optimal  feedback  controllers  for  actual  HVAC 
systems.  The  equations  which  accurately  describe  the 
controlled  processes  are  often  nonlinear  and  one  approach 
has  been  to  linearize  them,  to  periodically  update  the 
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parameters  of  the  linearized  model  by  relinearizing  the 
nonlinear  terms[4]  or  by  using  a  recursive  least  squares 
algorithm  to  identify  parameters  from  data[2],  and  then  to 
solve  the  matrix  Riccati  equation  for  near-optimal  control 
in  the  next  time  interval.  This  approach  requires  rather 
sophisticated  computations  at  each  update  and  is  inherently 
limited  in  its  optimization  look-ahead  time  because  of  the 
limited  duration  of  accuracy  for  the  linearized  model. 

We  propose  that  a  better  approach  is  to  explicitly 
recognize  that  many  of  the  HVAC  processes  can  be  accurately 
modelled  as  bilinear  processes,  where  net  heat  transfer  rate 
is  related  to  the  product  of  pumping  rate  and  temperature 
difference.  This  is  true  for  solar  collectors,  heat  pumps, 
and  building  heat  delivery.  Thus,  if  optimal  control  can  be 
computed  on  a  basis  of  a  bilinear  model  that  remain  globally 
valid,  the  need  for  frequent  model  updating  and 
reoptimization  is  eliminated  and  look-ahead  times  are  not 
limited  by  model  validity. 

The  solutions  of  bilinear  control-optimization  problems 
are  easily  obtained  for  simple  (first  order)  systems.  For  a 
solar  heating  system  with  one  storage  tank,  Auslander  et 
al[5]  have  shown  that  simple  energy  considerations  dictate 
that  the  building-temperature  control  problem  can  be 
decoupled  from  the  collector  control  problem.  Any  good 
regulator  that  holds  building  temperature  close  to  its 
minimum  specified  level  is  nearly  optimal.  Moreover,  the 
optimal  collector  pump  controller  is  bang/bang  (or  on/off), 


when  collected  solar  energy  is  to  be  maximized  and  parasitic 
(pumping)  losses  are  not  subtracted.  This  follows  from  the 
fact  that  the  Hamiltonian  is  linear  in  control  variable 
(which  is  the  collector  heat  removal  factor,  a  monotonic 
function  of  pumping  rate)  and  singular  arcs  can  be  shown  not 
to  exist. 

With  the  same  system,  but  with  parasitic  losses 
considered,  Winn  and  Hull  [1]  have  shown  that  the 
Hamiltonian  becomes  convex  and  pumping  rate  becomes  a 
continuous  function  of  tank  and  collector  outlet 
temperatures.  There  is  some  doubt  that  pumping  power  should 
be  considered  as  a  loss,  however,  since  with  proper 
placement  of  components  this  power  will  also  be  a  heat  input 
to  the  system,  hence  will  offset  a  cor  responding  amount  of 
auxiliary  heat  input. 

A  more  practical  solar  heating  system  for  northern 
lattitudes  will  probably  turn  out  to  be  the  series 
solar-assisted  heat  pump  system,  however,  because  of  the 
neccessar i ly  large  difference  between  the  building- loop 
inlet  and  ambient  temperatures.  If  two  storage  tanks  are 
included,  as  shown  in  Fig.  2.1,  the  system  offers  the 
further  attractive  option  of  efficient  day-time  solar 
collection  with  a  relatively  cool  collector  feeding  a  low7 
temperature  tank  and  efficient  night-time  (offpeak)  heat 
pumping,  with  a  good  coefficient  of  performance,  to  the 
high-temperature  tank. 
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The  optimal  control  problem  now  becomes  much  more 
difficult,  however,  because  while  the  collector  and  building 
controls  are  still  decoupled,  the  optimal  heat  pump  control 
sequence  must  depend  on  both  tank  temperatures  and  on  the 
anticipated  building  heating  load  and  solar  intensity 
patterns . 

The  formal  optimization  is  difficult,  since,  while  the 
Hamiltonian  is  still  linear  in  the  controls,  it  can  be  shown 
that  singular  arcs  are  possible  with  this  system  and, 
indeed,  will  be  involved  in  the  optimal  control  sequence. 
Gunewardana  et  al  [6]  have  successfully  used  a  numerical 
search  method  to  compute  optimal  preview  control  sequences 
for  a  similar  nonlinear  system. 

In  the  following  section  we  obtain  the  optimal  control 
sequences  by  numerical  search  method,  we  then  postulate  that 
the  near-optimal  controller  could  feasibly  be  implemented 
with  a  microprocessor  that  (1)  preselects  a  set  point  that 
depends  in  a  relatively  simple  way  on  the  weather  forecast 
(2)uses  an  adaptive  controller  of  which  the  current  value  is 
determined  from  the  previous  one  (more  will  be  said  about 
this  controller  in  the  next  chapter) 


2.2  The  system  model. 
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2.2.1  The  heating  load. 

Using  the  example  13-2  of  [10],  the  house  is  assumed  to 
be  constructed  as  follows: 

Exterior  walls: 

4"  common  brick 
1/2  "  plywood 
2x4  studs 
R- 1 1  insulat ion 
1/2  "  plasterboard 

Floor  construction  over  vented  crawlspace 
25/32  "  hardwood  finish  flooring 
building  paper 
1"  plywood  sub-floor 
airspace 

R- 1 1  insulation  (applied  to  underside  of  joists) 
Windows:  storm  windows  (4) 

Exterior:  1-1/2  "  solid  core  door 

Ceiling  construction  with  vented  attic  space  above 
1/2"  plasterboard 
R- 1 9  insulation 

The  house  is  51'x27',  and  has  one  wood  exterior  door 
3'x6'-8",  and  one  double  glass  wood  frame  sliding  patio  door 
with  1/4"  air  space,  7'x6'-8'.  The  heat  load  calculation  can 
be  summarized  as  follows: 
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u 

( Btu/hr  f  t 2  °F ) 

A 

(  f  t !  ) 

ub 

(Btu/hr  °F ) 

Exterior  walls  0.07 

974 

68.18 

Windows  and  sliding  patio  doors 

double  0.65 

47 

30 . 55 

single  0.56 

207 

115.92 

Exterior 

slab  doors  0.49 

20 

9.8 

Floors  0.07 

1377 

96.39 

Ceilings  0.05 

1377 

68 .85 

Total  389.69 

For  simplicity,  infiltration  loss  is  ignored. 

If  we  designate  a^  as  the  overall  heat-loss  coefficient  of 

b 

the  building,  we  have: 

a  =389.69  Btu/hr °F=740 . 06  KJ/hr°C  (2.2-1) 


2.2.2  Collector  model. 
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The  collector  assumed  is  a  flat  plate  collector.  The 

rate  of  solar  energy  collection  Q  is  given,  as  previously, 

c 

by  equation  (1.4-3) 

Qc=AcF  (u3  )  [HTar-UL(Ts-Ta  )  ] 

Let  S=H^ar,  we  have: 

Qc=AcF(u3 ) [S-U  (Ts-Ta) ]  (2.2-2) 

Using  the  numbers  for  the  example  collector  given  in  [1],  we 
have 


Ac= 1 00  m 
ar  =  Q  .  84 


F  =0.867 


Ul=12.  KJ/(hr*m2*°C) 

The  collector  has  to  be  tilted  to  optimize  the  reception  of 
solar  radiation,  but  only  radiation  data  incident  on  a 
horizontal  surface  are  available.  Hence,  the  following 
parameters  are  needed  to  process  these  data  in  order  to 
obtain  values  for  Km  (see  Appendix  I). 


The  collector  is  placed  due  south,  Hence  7=0 


o 


The  location  is  in  Edmonton,  of  which  lattitude  is  54 


o 


The  tilt  slope  of  the  collector  is  64°(=lattitude+10u) 


,o 


2.2.3  The  storage  tank  model. [13]. 

There  are  two  ways  of  storing  heat,  either  as  sensible 
heat  or  as  latent  heat. 
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Energy  may  be  stored  as  sensible  heat  in  liquid  (e.g. 
water)  or  solid  medium  (e.g.  rock)  while  selected  salts  are 
used  to  to  store  latent  heat  in  the  form  of  heat  of 
crystallization.  It  is  sensible  heat  that  has  been  widely 
and  reliably  used,  although  a  lot  of  research  is  being  done 
to  understand  more  about  the  use  of  latent  heat.  The  amount 
of  energy  that  can  be  stored  as  sensible  heat  in  a  storage 
tank  is: 

Q  =mC„AT  (2.2-3) 

t  P 

» 

where, Q=heat  capacity  of  system. 
m=mass  of  storage  medium. 

Cp-Specific  heat  of  storage  medium. 

AT=Tempe ra t. ure  range  between  which  the  tank  operates. 
Losses  of  sensible  heat  stored  in  a  tank  are  due  to 
conduction  through  the  temperature  difference  between  the 
tank  and  the  surrounding  environment. 

6  =a  ( T , -T 2 )  (  2 . 2-4  ) 

1  t 

Where,  aris  the  heat  loss  coefficient  of  the  tank 
T,  is  the  tank  temperature 
T2  is  the  environment  temperature 
There  are  two  storage  tanks  in  our  system,  one  of  which 
operates  at  low  temperature,  and  is  called  the 
low- temperature  tank  or  low-side  tank.  The  other  tank,  which 
works  at  higher  temperature,  is  called  the  hi gh- temperature 
tank  or  high-side  tank. 

The  low-side  tank  is  used  to  collect  the  available 
solar  energy  from  the  collector.  It  is  often  operational  at 
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a  higher  temperature  than  the  ambient.  By  absorbing  all 
available  energy  of  the  collector  with  only  small  increase 
in  its  average  temperature,  the  low-side  tank  helps  keep  the 
average  temperature  of  the  collector  at  lower  value  than 
would  be  the  case  if  the  collector  were  connected  directly 
to  hi gh-tempera ture  storage  tank,  thus  increasing  the 
collector  efficiency.  Moreover,  the  temperature  of  the 
low-side  tank,  which  is  higher  than  the  ambient  temperature 
because  of  the  stored  solar  energy,  also  helps  boost  the 
coefficient  of  performance  of  the  heat  pump. 

The  energy  balance  at  rhe  low-side  tank  is  : 


C  T  =  Q  -  Q  -  Q, 
c  c  vc  t  hp 

where  C  =low-side  tank  heat  capacity, 


(2.2- 


Qc i s  the  rate  at  which  energy  is  drawn  from  the 


collector.  According  to  equation  (2.2- 


) 


Qc=AcF(u3)[S-UL(Tc-Ta)]  (2.2-6) 

here  T  is  replaced  by  Tr ,  temperature  of  the  low-side  tank. 


6  is  the  heat  loss  of  the  tank 
t 


(2.2-7) 


V£c(VV 

Q,  is  the  rate  at  which  heat  is  taken  out  of  the 
xhp 

low-temperature  tank  by  the  heat  pump.  An  expression  for  Q 
will  be  found  in  the  next  section. 


hp 


The  energy  balance  at  the  high-side  tank  is: 


t  *  f 


chV  %  V  Qb 

where  C^=high-side  tank  heat  capacity, 

o'  is  the  heat  loss  of  the  tank 

t 

Q  '  =a,(  T,~T  ) 

t  h  h  a 


(2.2-8) 


(2.2-9) 
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Q  is  the  heat  delivered  to  the  building 

L/ 

Q^u,  (Th"Tb)  (2.2-10) 

where  £  =heat  exchanger  coefficient 

T 

Q.  is  the  heat-pump  output  (see  next  section). 

hp 

According  to  Lof  and  Tybout  [15],  the  optimum  heat  storage 
is  in  the  range  of  50  to  75  Kg  of  water  per  one  square  meter 
of  the  collector  area.  The  value  of  75  Kg/m2  is  used  here. 


2.2.4  Heat  pump. 

Using  a  commercially  available  3-ton  heat  pump[14], 

which  is  manufactured  by  the  Carrier  Corporation,  a  graph  of 

the  COP  with  respect  to  the  temperature  difference  between 

the  evaporator  and  the  condenser  is  plotted  in  Fig. 2. 2. 

The  coefficient  of  performance  can  be  approximated  by  a 

straight  line  whose  equation  is: 

COP(Th,Tc)  =  1  +  (COPmax-1)(1-(Th-Tc)/Tmax)  (2.2-11) 

Where,  COP-,  =3.5  and  T  =45.°C 
max  max 

If  we  designate  u2  as  the  heat-pump  electrical  input,  the 
rate  at  which  heat  is  taken  out  of  the  low-side  tank  by  the 
heat  pump  is: 


Q,  =u2  [COP(T,  ,  T  )  -  1  .  ]  . 
ho  he 

The  heat  pump  output  is  calculated  as  follows: 


(2.2-12) 


VU2C0P<Th'V- 


(2.2-13) 
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2.2.5  Circulation  pumps. 

According  to  experiments,  the  collector  flow  rate 
should  be  in  the  neighborhood  of  60  liters/hr  per  one  square 
meter  of  collector  area  [12].  Actually,  the  value  of  u3, 
collector  flow  rate,  is  not  needed  as  far  as  equation 
(2.2-2)  is  concerned. 

The  building  loop  pump  capacity  is  calculated  as 
follows: 

The  heat  load  of  the  building  is: 

Q  =  a  ,  ( T,  -  T  )  KJ/hr  (2.2-14) 

1  b  b  a 

where  T  =building  temperature=20  C, 

T^  =  average  ambient  temperature. 

Take  T  =-12°C  [I2](average  of  January) 

3 

The  maximum  heat  delivered  by  the  high-side  tank  is: 

2b=i;u’n>ax(VV  KJ/hr 

where  £=0.8 


T,  =average  temperature  of  the  hich-side  tank. 

h 

Assume  that  T^  is  approximately  equal  to  24°C. 

In  order  for  the  building  heat  load  ro  be  fully  met  by  the 
heat  from  the  high-temperature  tank  under  these  conditions, 
the  building  loop  pump  capacity  must  be  equal  to: 
“■»ax'ab(Tb-Ta)/t(Th-Tb) 

=740. 06x32. /0.8x(24. -20. )=7400.  Kg/hr 


2.3  The  Steady-State  Optimal  Temperatures 
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We  assume  that  an  effective  load  regulator  holds  the 
building  temperature  constant  at  a  specified  T^using 
auxiliary  heat  if  necessary.  The  collector  pump  rate 
controller  is  a  bang  bang  optimal  controller: 


u3=U3maxt 1 .+sgn(S-UT (T  -T  ))]/2. 

Li  c  a 


I  .  e  .  , 


(2 . 3- la  ) 


Ta» 


(2.3-1b) 


If  S  T^ ) > 0 .  u3-u3max 

If  Q  =S-U(T  -T  )<0.  u  3  =  0  . 
o  L  C  a 

The  dynamic  equation  for  the  two  tanks  are  then: 

CUT  =-  £  u ,  ( T,  -T,  )  -  a  ( T  -T  )  +u  2  COP  ( T  ,  T  ) 

hh  hb  hha  h  e 

C  T  =~u  2 [ COP ( T  , T  )- 1 ]-a„ (T  -T a)+A  F(u3)[S~Ut{T  - 

C  C  hC  CCdC  LC 


where  the  Coefficient  of  Performance  is: 

COP  ( T  , T  )  =  1 .  +  ( COP  -  1 . ) [ 1  - ( T  -T  )/T  3 

h  c  max  h  c  max 

and  where: 

T,  ,  T  ,T  ,T  are  hot  tank,  cool  tank,  ambient,  end  buildinc 

n  c  a  5 

temperatures . 

C , , C  are  tank  heat  capacities 

h  c 

u,,u2,u3  are  building  loop  pump,  heat  pump, and  collector 
pump  rates. 

£,a,  ,a  are  heat  transfer  coefficients 

h  c 

Ar,F  ,U  are  collector  area,  heat  removal  factor,  and  loss 
coefficient. 

The  required  auxiliary  heat  input  rate  to  the  building 
i  s 


£>  =aw(Tu"TJ-£ui  (Tu"TJ^0 

aux  b  b  a  h  b 

where  the  equality  holds  if  u,£  u,maxcan  satisfy  the 


equation,  i.e.  no  auxiliary  heat  is  used  if  the  building 
load  can  be  supplied  from  storage. 

The  steady-state  optimal  problem  can  be  stated  as  follows: 

Minimize  J=a,(T  ~T  )-£u,(T  -T,  )  +  u2 

b  b  a  n  b 

(J  represents  the  total  of  auxiliary  energy  and  the  input 
power  to  the  heat  pump) . 
subject  to: 


-£u, (T-T,  )-aL(T-T  )+u2COP(T  fT  )=0 

n  b  n  n  a  h  c 


(2.3-1  ) 

-u  2 [ COP (T  , T  )  -  1  ]  -  a  (T  -T  ) +A  F ( u  3 ) [S-U  (T  -T  ) ]  =  0  (2.3-2) 

n  C  CCdC  LCa 

(2.3-3) 


-a  (T  -T  )+  u,  (T,  -T,  )<0 

b  b  a  h  b  v 

u  1  ~u  ’max 


<0 


u  2  u 2 max  . 

Establishing  Lagrange  function  and  then  using  the 
KUHN-TUCKER  conditions  to  solve  this  problem[9], 
the  steady-state  optimal  high-side  temperature  will  be 
proved  to  be: 


(2.3-4) 

(2.3-5) 


h  b  b  a 

°,u°  are  found  from  (2.3-1)  and  (2.3-2) 


)  +T 

max  ’  1  b 


Proof : le t  the  Lagrange  function  be  defined  as  follows: 
L=J+p,*lhs(2.3-1)+p2*lhs(2.3-2)+p3*lhs(2.3-3) 


+p4*lhs(2.3-4)+p5*lh'S(2.3-5) 


where  lhs (2.3—1 )  is  the  left  hand  side  of  the  equation 
(2.3-1) 

lhs (2.3-2)  is  the  left  hand  side  of  the  equation 
(2.3-2) , etc  .  .  . 


The  KUHN-TUCKER  conditions  are: 
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c'L/^T,^  =  0  ,  dL/c)Tc  =  0  , 

^L/^u, =0 , dL/^u  2  =  0, 

Constraint  equations  (2.3-1),  (2.3-2),  (2.3-3), 
(2.3-4),  (2.3-5), 

p. *lhs ( 2 . 3-i ) =0 . ,  i=1,2,...,5,  and 
P.^0,  i  =  3, 4 , 5 


Applying  the  KUHN-TUCKER  theorem,  we  obtain  the  following 
condi t i ons : 


1 )  =  - {u,+p, [-{u,-a^-u2* (COPmax  -  1 ) /Tmax) 

+  Pj  [  u2  *  (COPma)(  -  1  )  /Tmax  ]+ps[(u,]  =  0.  (2.3-6) 

2) SL/aTc  =  p,[uJ*(COPmax-1)/Tmax]+p1*(-u2*(COPIiax-1)/Tmax 

-a. -  A  F  ( u  2  )  U  ]  =  0  (2 .3-7 i 

w  c  L 

3  )  <5L/<5u  2  =  1  +p  i  COP  (  T  ,  T  )  -p  2  *  [  COP  (  T  ,T  )-1  ]+p5  =  0  (2.3-8) 

n  c  h  u 

4  )  c>L/^u  ,  =  -  HT  -T.  )+Pl  [-HT.  -T,  )  ]+p3  [  £(T  -T  )  ]+p4  (2.3-9) 

n  D  n  b  h  b 

=0.  (2.3-9) 


5) The  constraint  equations 

(2.3-4)  ,  (2.3-5)  . 

6)  p1*lhs(2.3-1)=0 

7)  p2  *lhs ( 2 . 3-2 )  =  0 

8)  p3*lhs(2.3-3)=0 

9)  p4 *lhs ( 2 . 3-4  )  =0 

10)  p 5  * lhs ( 2 . 3- 5 ) =0 

11) no  restrictions  on  p, 

12) no  restrictions  on  p2 

1 3) p3>0 


(2.3-1  )  ,  (2.3-2)  ,  (2.3-3)  , 

(2.3-10) 

(2.3-11) 

(2.3-12) 

(2.3-13) 

(2.3-14) 

(2.3-15) 

(2.3-16) 

(2.3-17) 
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1 4 )p4>0 

(2.3-18) 

1  5  )  p  5  ^0 

(2.3-19) 

Equation  (2.3-12)  gives: 

p3*[-a,  (T.  -T  )  +  $u,  (T  -T  ) ]-0. 

D  D  a  h  b 

(2.3-12) 

Case  1  :  p3  =0 . 

(2.3-20) 

Equation  (2.3-13)  gives:  p 4 [ u , -u , mdx ] =0 . 

(2.3-13) 

Case  1.1:  p4=0 . 

(2.3-21  ) 

Equation  (2.3-14)  gives:  p5 [ u2 -u2max  ] =0 . 

(2.3-14) 

Case  1.1.1:  p5  =  0 

(2.3-22) 

From  (  2 . 3-20  )  , ( 2 . 3-2 1 ) , ( 2 . 3-9 ) , we  obtain: 

♦o 

it 

i 

• 

(2.3-23) 

Substituting  (2.3-23)  into  (2.3-6)  and  (2.3-7),  we  find  that 
there  are  no  values  of  p2  that  can  satisfy  both  (2.3-6)  and 
(2.3-7)  . 


Therefore,  this  case  is  excluded. 

Case  1.1.2:  u2=u2(nax  (2.3-25) 

This  case  is  excluded  for  the  same  reason  as  the  above. 

Case  1.2:  u,=u,max  (2.3-27) 

Case  1.2.1:  ps=0.  (2.3-28) 

From  (2.3-8): 

P,  =  (p2*[COP(Tu ,Tr  ) -  1 ] -  1 ) / COP ( T  , T  ) .  (2.3-29) 

Replacing  (2.3-29)  into  (2.3-7)  ,and  rear  rang  i  ng , we  obtain: 
p.-[u1.(C0Pmax-l’)/(TmaxC0P(Th,Tc))]/[-uI*(C0pmax-1) 

/(T  C°P(T  T  ))-a  -A  F(us)u  ].  (2.3-30) 

max  he  cc  |_ 

Equation  (2.3-6)  is  violated  when  we  substitute  the  values 
of  p,  ,  p2 ,  p 3  which  are  described  respectively  by  equations 


(2.3-29),  (2.3-30),  (2.3-20). 

Therefore,  this  case  is  also  excluded. 
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Case  1.2.2:  u2=u2max.  (2.3-31) 

Rewriting  (2.3-7),  we  obtain: 

p,-p,  =  p![ac+A(;F(uJ)UL]/  [u,max*<COPmax-1)/Tmax]  (2.3-32) 

Because  OOPmax>0  and  p2^0,  Hence  from  (2.3-32) 

p,-p2^0.  (2.3-33) 

Computing  p5  from  (2. 3-8), we  have 

Ps  =  - 1 -p , COP ( T  , T  ) +p2 [COP(T  ,T  )  -1]. 

n  C  h  c 

Or, 

Ps=- 1-p2-COP(T  , T  ) [ p , -p  2 ]  (2.3-34) 

n  c 

From  (2.3-33),  (2.3-34)  and  for  COP(T  ,T  )>0  ,  we  derive: 

n  c 

ps<0 

The  above  condition  violates  the  constraint  inequality 
(2.3-19).  Hence,  this  case  is  excluded. 


Case 

2. :  a  ( T  -T  )-£u, (T  -T  )=0 

b  b  a  h  b 

(2.3-35) 

Case  2.1:  p 4  =  0 

(2 . 3-36) 

Case  2.1.1:  p5  =  0 

(2.3-37) 

From 

equations  (2.3-9)  and  (2.3-36),  we  obtain: 

p3= 1 +p, 

(2.3-38) 

Repla 

cing  (2.3-38)  into  (2.3-6),  we  get: 

p,  [-a 

^-u2  *  (C0Pmax~ 1  )/Tmax  [u2  *  (c°Pmax  - 1  )/Tmax  =  0 

(2.3-39) 

From 

(2.3-7)  and  (2.3-39),  we  obtain: 

p, =p2=0 

(2.3-40) 

Replace  (2.3-37)  and  (2.3-40)  into  (2.3-8) 

1=0  (Contradiction). 


This  case  is  excluded. 
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Case  2.1.2:  u2=u2mfl  Y 

(2.3-41 ) 

This 

case 

is  also  excluded,  the  proof  is  exactly 

like  that 

of  case  1 

.2.2. 

Case 

2.2:  u,=u,max 

(2.3-42) 

Case  2.2.1:  u2  =u 2ma  v 

(2.3-43) 

This 

case 

is  excluded  for  the  same 

reason  as  that 

of  case 

2.1. 

2  or 

1.2.2. 

Case  2.2.2:  p5  =  0 

(2.3-44) 

From 

(2.3 

-42)  and  (2.3-35),  as  was 

claimed 

T?  =  aJTu"ri;>  )/U*u 

n  b  b  a 

’max ■ * 

(2.3-45) 

T°,u°  are  found  from  (2.3-1),  and  (2.3-2)  by  solving  a 
system  of  equations.  The  analysis  to  find  existence 
conditions  for  T° ,  T°  is  omitted  because  the  computations 
are  very  involved.  Moreover,  the  suboptimal  controller  is 
not  smart  enough  to  take  these  conditions  into  account.  This 
is  one  area  that  needs  further  studies.  However,  we  note 
that  when  S  equals  zero  the  steady- state  optimal 
temperatures  for  both  tanks  do  not  exist. 

The  proof  is  as  follows: 

since  T°>Uh ( because  of  stored  solar  energy), 

C  a 

hence,  S-U^ ( T°-Td ) <0 . 

This  implies  F(u3)=0. 

From  (2.3-2),  we  obtain: 


since, 
hence , 


u°2  [COP(T°,T°)-1  .  ]=-ac(T°-Ta  )  . 

COP(T° , T° ) -  1 . >  0 .  , 

h  c 


Because  u2  is  the  heat  pump  electrical  input,  it  cannot  take 


36 


on  negative  values. 

Hence,  u°  =  0.  (2.3-46) 

From  (2.3-46),  (2.3-42),  (2.3-45)  and  (2.3-1),  we  obtain: 


Th  ^ahTa+stUl'liaxT^//^a^+^u'max^  • 

There  are  no  values  of  T°  that  can  satisfy  both  the 

equations  (2.3-45)  and  (2.3-47). 

This  completes  the  proof. 


(2.3-47) 


2.3  Optimal  control  sequence. 

Although  the  steady- state  opt imal  control  is  of 
interest  in  its  own  right,  we  are  here  more  interested  in 
the  nature  of  the  heat  pump  control  sequences  that  are 
optimal  under  other  than  steady-state  conditions.  This  is 
because  the  solar  intensity  S  and  ambient  temperature  T  are 
never  fixed  for  a  long  time  although  they  are  expected  to 
change  gradually.  Therefore,  the  heat  pump  control  will  keep 
changing  to  move  the  system  operating  point  toward  the  new 
steady-state  operating  point  when  S  and  T  change.  The 
optimal  values  of  u2  calculated  as  above  will  not 
necessarily  be  optimal  under  these  conditions. 

The  steepest  descent  method  is  used  to  search  for  the 
optimal  values  of  the  heat  pump  control  sequence  that  will 
minimize  a  cost  index  which  is  a  function  of  auxiliary  heat 
and  electrical  energy  input  to  the  heat  pump. 
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2.3.1  The  steepest  descent  method[8]. 

The  gradient  of  a  function  is  the  vector  of  partial 
derivatives  of  the  function,  f,  with  respect  to  each  of  its 
variables.  The  gradient  has  a  very  important  property,  that 
is  the  function  value  increases  at  the  fastest  rate  when  we 
move  along  the  the  gradient  direction.  Therefore,  the 
gradient  direction  is  called  the  steepest  ascent  direction. 
Similarly,  the  negative  of  the  gradient  direction  indicates 
the  direction  of  steepest  descent. 

The  steepest  descent  method  is  the  method  for 
minimizing  a  cost  function  that  makes  use  of  the  direction 
of  the  steepest  descent.  In  this  method,  we  start  from  an 
initial  guess  point  X,  for  the  argument  of  the  function  and 
iteratively  move  towards  the  optimum  point  according  to  the 
formula : 


X .  =X .+  X. S  . 

-i  +  l  -i  i-i 

where  S_.is  the  negative  of  the  gradient  of  the  function  f, 
X.is  the  optimal  step  size  which  is  obtained  by  making  an 
one-dimensional  search  along  direction  S  . 

Because  the  gradient  magnitude  gets  smaller  as  it  is 
approaching  the  optimum,  the  rate  of  convergence  becomes 
smaller  and  smaller.  In  order  to  improve  the  convergence 
characteristics  of  the  steepest  descent,  the  search 
direction  needs  to  be  modified.  One  of  the  modified  steepest 
descent  methods  is  the  conjugate  gradient  method  which 
determines  the  current  search  direction,  based  upon  the 
previous  direction,  the  previous  gradient,  and  the  current 
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gradient  of  the  function.  This  method  has  been  shown  to 
converge  quadrat ically ,  and  it  will  be  used  to  determine  the 
optimal  control  sequence  of  the  heat  pump. 

Using  the  algorithm  suggested  by  Fletcher  and  Reeves  to 
minimize  general  functions  as  follows[8]: 

1)  Start  with  an  arbitrary  point  X, 

2)  Set  the  first  search  direction 


S, =-Vf (X, )=-Vf , 

3)  Find  the  point  X2  according  to  the  relation 

x2=x,  +  x1s1 

where  X^is  the  optimal  step  length  in  the  direction  S, 
Set  i=2  and  go  to  the  next  step. 

4)  Find  Vf..=Vf(X.)  and  set 


S^-Vf.+t  |V£.  I  VIVf.^1  MS 


1 


By  trial  and  error  process,  the  factor  multiplying  S/  can  be 
replaced  by  a  constant. 

5)  Compute  the  optimum  step  length  X  ^  in  the  direction  S , 


and  find  the  new  point: 


-i+1  -i+Xi"i 


6)  Test  for  the  optimality  of  the  point  X ..  . 

If  X  is  optimum,  stop  the  process.  Otherwise,  set  the 
value  of  i  equal  to  i+1,  and  go  back  to  step  4). 

The  program  used  to  find  the  optimal  path  of  this 
system  controller  is  listed  in  appendix  II  (the  weather 
conditions  used  are  hypothetical).  Specifically,  we  can 
divide  this  program  into  three  parts.  The  first  part  of  the 
program  is  to  establish  the  initial  trajectory  of  the 
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controller.  In  the  second  part,  the  gradient  of  the 
objective  function  is  evaluated,  and  then  an  one-dimensioned 
search  is  performed  along  the  direction  determined  by  the 
gradient.  The  search  is  stopped  in  part  three  of  the 
program . 

The  initial  control  trajectory  is  developed  as  follows: 
The  building  loop  pump  controller  u,  is  turned  on  to  such  a 
value  that  the  building  heat  load  is  fully  met  by  the  heat 
from  the  high-side  tank.  In  the  case  the  building  load 
cannot  completely  be  statisfied  even  with  the  maximum  value 
of  u, ,  the  rest  of  the  load  is  supplied  from  the  auxiliary 
heater.  The  heat  pump  controller  u2  is  set  to  such  a  value 
that  the  heat  removed  from  the  high-side  tank  will  be 
replenished  by  the  same  amount  of  heat  from  the  low-side 
tank.  If  this  cannot  be  realized,  the  heat  pump  is  turned  on 
to  its  maximum  value.  The  collector  pump  controller  is  of 
on-off  nature.  It  is  either  switched  on  to  the  maximum  value 
or  switched  off  completely,  depending  on  whether  or  not  the 
solar  energy  collected  is  more  than  the  heat  loss  from  the 
collector  to  the  ambient. 

Noting  that  the  heat  pump  controller  sequence  has  50 
sampled  values,  each  of  which  corresponds  to  a  time 
interval.  In  the  second  part  of  the  program,  the  gradients 
are  found  by  imagining  that  we  have  a  space  of  35 
dimensions,  each  represents  one  out  of  the  first  35  time 
intervals  of  the  controller  time  sequence.  The  gradient 
vector  can  be  thought  of  as  a  35-dimensioned  vector  whose 


40 


each  component  is  calculated  by  varying  the  corresponding 
value  of  u2  along  that  dimension.  The  last  15  sampled  values 
of  u2  are  used  to  bring  the  system  states  back  to  their 
initial  values,  thus  preventing  the  stealth  of  energy  from 
the  final  states  to  save  on  purchased  energy  (the  objective 
function ) . 

The  new  search  direction  is  calculated  from  the  current 
gradient  vector  and  the  old  search  direction.  By  experiment, 
the  coefficients  of  the  gradient  vector  and  the  old  search 
direction  are  chosen  as  -0.5  and  0.5,  respectively.  Here, 
the  so-called  one-dimen sioned  search  actually  means 
searching  along  one  direction,  which  is  determined  as  above. 

The  third  part  establishes  the  criterion  to  stop  the 
search . 

The  optimum  control  sequence  of  the  heat  pump  obtained 
from  running  the  above  program  is  plotted  in  Fig.2.3[7]. 


2.3.2  The  proposed  suboptimal  controller. 

Using  the  model,  which  is  developed  previously,  the 
numerical  optimal  path  of  the  heat  pump  controller  sequence 
was  found  by  the  steepest  descent  method.  However,  for 
practical  reasons,  it  is  very  difficult  to  implement  the 
optimal  controller  sequence.  Therefore,  an  implemen table 
suboptimal  controller  needs  to  be  developed. 


Various  avenues  have  been  explored  to  develop  a 
suitable  suboptimal  controller  for  this  system.  One 
promising  approach,  which  is  explained  in  detail  in  Appendix 
V  can  be  be  applied  to  a  general  bilinear  system.  However, 
more  work  needs  to  be  done  before  this  suboptimal  controller 
can  be  implemented  satisfactorily. 

We  have  also  devoted  much  time  to  another  algorithm. 
That  is,  using  the  constant  COP  path  and  the  discrimination 

curve  in  the  state  plane  of  T  and  T  (Fig.  2.4),  the  system 

h  c 

can  be  moved  from  its  current  operating  point  to  the 
corresponding  steady-state  optimum  point  when  the  weather 
changes.  For  example,  the  current  weather  is  fairly  good  and 
the  system  operating  point  is  at  point  A  in  the  state  plane 
(Fig.  2.4).  Assuming  that  the  predicted  weather  is  going  to 
be  extremely  bad,  thus,  based  upon  this  prediction,  a  new 
steady-state  optimal  point  B  (Fig.  2.4)  is  determined.  The 
problem  is  how  and  when  to  get  to  the  new  point  from  the  old 
one.  One  proposed  way  is  to  adjust  the  heat  pump  such  that 
its  COP  (coefficient  of  performance)  is  constant,  thus 
moving  the  system  along  a  constant  COP  curve  in  the  state 
plane.  The  reason  for  this  is  to  keep  the  heat  pump 
operating  at  relatively  good  values  of  COP  under  adverse 
weather  conditions  while  waiting  for  the  good  weather  to 
come.  When  the  discrimination  function  changes  sign,  heat 
pump  will  be  turned  to  its  maximum  value  to  reach  the  new 
optimum  point  (discrimination  curve  is  a  curve  which  passes 
through  the  new  optimum  point,  with  heat  pump  set  to  its 
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maximum).  The  timing  of  the  heat  pump  is  adjusted  such  that 
the  new  optimum  point  is  reached  just  when  the  good  weather 
conditions  arrive.  But  there  are  doubts  that  such  a 
trajectory  always  exists  for  various  weather  conditions,  let 
alone  being  near-optimal.  It's  also  impractical  to  expect 
such  an  abrupt  change  of  the  weather.  Moreover,  the 
computation  will  be  very  much  involved,  therefore  destroying 
any  prospect  of  its  implementation  on  a  microprocessor . 

Based  upon  the  results  found  from  running  the  above 
program  with  real  conditions  of  weather,  and  upon  the 
results  of  a  process  of  trial  and  error,  we  propose  a 
suboptimal  preview  adaptive  controller  sequence  as  follows: 

u 2 ( n+ 1 )  =  0 . 5*u 2 ( n ) +0 . 5*K ( n ) * ( 1 . +  T,  -T,  (n) )  (2. 4. 2-1) 

ns  h 

where , 

u2(n+1)=the  next-step  value  of  the  controller, 
u2(n)=the  curent  value  of  the  controller, 

K(n)=an  adaptive  gain  changed  every  step 

=ratio  of  the  total  high-side’s  heat  loss  to  the 
heat  pump  Coeffient  of  Performance  at  step  n, 

=[{u, (n ) (T  (n)-T  )+a(T(n)-T  )]/COP(T  (n),T  (n)) 

h  b  h  h  a  h  c 

T  =The  steady-state  optimal  temperature  of  the 

hs 

high-side  tank  (calculated  by  using  equation  (2.3-45)  with  TQ 
being  the  predicted  average  daily  ambient  temperature). 

T^(n)=The  current  temperature  of  the  high-side  tank. 

In  the  following,  we  will  prove  that  the  above 
suboptimal  controller  tends  to  bring  the  high-side 
temperatures  back  to  its  steady-state  optimal  values. 
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Assuming  that  T  (n)<T  (2. 4. 2-2) 

h  hs 

u  2 ( n ) >K ( n )  (2. 4. 2-3) 

According  to  definition,  K(n)  is  meant  to  be  the  required 
heat  pump  electrical  input  to  keep  the  high-side  temperature 
unchanged  (see  equation  (2.3— la)  )  . 

From  (2. 4. 2-3),  we  obtain: 

u2 (n)>K(n)=[£u, (n) (T(n)-T)+a(T(n)-T  ) ]/COP(T(n) ,T  (n) ) . 

h  b  h  h  a  h  c 

From  (2.3-1a)  and  the  above  condition,  we  derive: 


Hence ,  T  ( n+ 1 ) >T  ( n ) 

h  h 

Therefore,  the  current  value  of  u2  tends  to  bring  the 
high-side  temperature  back  to  its  steady-state  optimal 
value  . 

On  the  other  hand,  if  T  (n)<T,  (2. 4. 2-4) 

h  hs 

and  u2(n)<K(n)  (2. 4. 2-5) 

From  (2. 4. 2-1),  (2. 4. 2-4),  and  (2. 4. 2-5),  we  obtain: 

u  2 ( n+ 1 ) >u2 ( n ) . 

Hence,  as  long  as  (2. 4. 2-4)  holds,  values  of  u2  will  be 
brought  up  until  u2  is  greater  than  K  (the  previous  case). 

If  (2. 4. 2-4)  is  violated  but  (2. 4. 2-5)  still  holds  we  go  to 
the  next  case. 


I  f 
and 


T  ( n ) >T 

h  hs 

u  2 ( n ) <K ( n ) 


(2. 4. 2-6) 
(2. 4. 2-7) 


(2. 4. 2-7)  implies  that 

T  ( n+ 1 ) <T  ( n ) . 

h  h 

The  current  value  of  u2  tends  to  stabilize  the  system  by 
bringing  the  high-side  temperature  back  to  its  steady-state 
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T  (  n  )  >T 

h  hs 


(2.4  .2-8) 
(2. 4. 2-9) 


optimal  temperature. 

For  the  last  case, 
if 

and  u  2 ( n ) >K ( n ) 

From  (2. 4. 2-8),  (2. 4. 2-9),  and  (2. 4. 2-1),  we  obtain: 

u  2  (  n+ 1 ) <u  2 ( n )  . 

Thus,  as  long  as  (2. 4. 2-8)  still  holds  values  of  u2  will  be 

decreased  until  u2  is  less  than  K  (return  to  the  previous 

case).  However,  if  (2. 4. 2-8)  is  violated  while  (2. 4. 2-9) 

still  holds,  return  to  the  first  case. 

In  conclusion,  regardless  of  the  current  values  of  T 

and  u2,  the  suboptimal  controller  tends  to  bring  the 

high-side  temperature  back  to  its  steady-state  optimal 

values.  The  reason  for  this  is  to  minimize  the  auxiliary 

heat  input  to  the  building. 

S  ince  ,  i f  T,  ( n ) =  T 

h  hs 


u  i  =  [  a  ,  (T,  -T,)]/U(Ths-Tb)]=uimax 


So , 


b  ’  b  a 


*ax<Ths-V  =  °. 


However,  if  more  information  on  T  is  taken  into  account  so 

c 

that  u2  is  minimized  at  the  same  time  with  £>,„,  there  can 
be  further  improvement  on  the  system  performance.  The 
approach  described  by  Appendix  V  may  be  useful  in  doing  so. 
The  stability  and  the  near-optimality  of  this  controller  are 
fully  explored  in  the  next  chapter. 
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Fig. 2.1:  The  System  Configuration 
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Fig.  2. 2:  Heat  Pump  Coefficient 
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Fig. 2. 3:  Optimum  Control  Sequence 
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Fig. 2. 4:  State  Plane 
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3.  SIMULATION  and  THE  SUB-OPTIMAL  CONTROLLER 


<; 


3.1  The  stability  of  the  suboptimal  controller. 

For  a  control  system,  stability  is  often  considered  as 
one  of  the  most  important  criterion  to  be  investigated. 
There  are  many  methods  such  as  the  Nyquist,  the 
Routh-Hurw i t z ,  the  Lyapunov  method, etc.,  that  can  be 
employed  to  analyse  the  stability  performance  of  a  system. 
In  the  following,  the  Lyapunov  method  will  be  used  to 
examine  the  stability  of  the  solar-assisted  heat  pump 
system . 

The  dynamic  equations  of  the  system  are  as  follows: 


x  =  A  x  +  (B  x  +  C)u  +(_D^  x  +  ID9)u^  +  _E  F(u^)  +  H 


where 


x. 

T  ~ 

1 

h 

x^ 

T 

2  _ 

c 

A  = 


B 


-(COP  -1) 
max 


max 


(3.1-1) 
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COP 

max 


1-COP 
_ max 

C 

c 


h  - 

? 

ch 

1 

o 

°2  - 

"5Tb  " 

Ch 

0 

0 

0 

_ 

E  = 


0 

.  rS-UT  (T  -T  )  , 
A  [  Leal 
c 


C 


Equation  (3.1-1)  has  the  form: 


x  =  f^(x,u1,u2,F(u3)) 


(3.1-2) 


In  the  neighborhood  of  the  steady-state  optimal  point,  define: 


6  x  =  x 


X  4- 

— opt 


(3.1-3) 
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6u 


1 


l,opt 


(3 . 1-4) 


6u 


2 


2  ,opt 


(3.1-5) 


6F(u3) 


F(u3)  - 


F(u,  ) 
j  ,opt 


(3.1-6) 


where  x  ,  u,  ,  u„  ,  u„  are  the  steady-state  optimal 

—opt’  1 ,  op  t  2, opt’  3, opt  J  r 

values  of  _x,  u^,  u^,  u^,  respectively. 

Substituting  (3.1-3),  (3.1-4),  (3.1-5),  (3.1-6)  into  (3.1-2), 
we  obtain: 


>x  =  f  (x  ^  +  6x  u1  ^  +  6u1  ,  u„  ^  +  6u  , 
— - opt  —  l,opt  1  ’  2, opt  2 


F(u,  )+6F(u  )) 
3, opt  i 


Linearizing  the  above  equation  around  the  optimum  steady  state 
point,  and  retaining  only  the  first-order  terms,  we  get: 


3f 


6  x  =  f  (  x 

u,  ,  ur 

F(u 

_  )) 

+  ~TTr 

ox  + 

—  —  —opt 

1  ,opt  / 

- , °pt ’ 

3 ,  opt 

9j* 

opt 

3f 

3f 

3f 

9u-^ 

6u^ 

+ 

3  2 

6u2 

+ 

"W(u3) 

6F(u3) 

(3.1 

opt 

opt 

opt 

where,  (...)  indicates  the  expression  (...)  is  evaluated  at  the 

opt 

optimum  point. 


From  (3.1-1),  we  have: 


3f 

3x 


=  A  +  _B  U2 


,opt 


opt 
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3u. 


opt 


D,  x  +  D 
— 1  —opt  —2 


3f_ 

^F(u3) 


=  B 

opt 


x  +  C 
—opt  — 


3f 

*3F(u3) 


=  E 


x 


opt  —opt 


Since,  f  (x  ,  u.  .  u  .,  F(u„  ))  =  x  =  0.,  (3.1-7)  can 
-  opt  !»opt-  2, opt  3, opt  -opt 


be  rewritten  as: 


<5 x  —  (A  +  B  u  )  6x  +  (D.  x  +  D,. 

—  —  —  2,  opt  —  —1  —opt  — ^ 


(B  x  +  C) 5u0 
- opt  —  2 


+  E_ i  6F (u3) 
I  x 

—opt 


(3.1-8) 


Here,  5F(u„)  can  in  principle  depend  on  x0  ^  =  T°,  since  F(u„)  is  a 
’  3  2 ,  opt  c  3 

bang-bang  controller  which  changes  its  value  when  S-U  (T  -T  )  =  0. 

J_»  C  3. 


If,  for  given  S  and  T  ,  T°  does  not  satisfy  the  above  switching 

3.  c 

condition,  i.e.  S-UT  (T°-T  )  4-  0,  there  will  be  a  neighborhood  of  T 

Lea  c 

values  which  do  not  satisfy  it  either.  Hence,  6f(u3)  =  0  in  the 

neighborhood  of  x 

°  —opt . 
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However,  if,  for  given  S  and  T  ,  T°  satisfies  the  switching 

3  C 

condition,  i.e.  S-UT (T°-T  )  =  0,  we  have: 

Lea 


x 

—opt 


0 

0 

A  [S-UT  (T°-T  )  ] 
c  Lea 

i 

o 

_  J 

C 

c 

=  0 


In  short,  e!  6F(u^)  is  always  zero, 


'  x 

—opt 


So,  (3.1-8)  becomes: 


6x  =  (A  +  B  u_  )  5x  +  (D-.  x  +  D_)  6u.,  + 
—  —  —  2,  opt  —  —1  —opt  —2  1 


(B  x  +  C)  <5u 
- opt  —  2 


(3.1-9) 


According  to  Appendix  IV,  the  homogeneous  system 


5x 


=  (A  +  B  u  )  fix 

-  -  2, opt  - 


(3.1-10) 


has  negative  eigenvalues  for  any  given  values  of  u9  q  .  Hence,  the 
system  described  by  equation  (3.1-10)  is  asymptotically  stable  [3]. 

System  (3.1-9)  reduces  to  system  (3.1-10)  when  6u^->  0.  Hence,  according 

2 

to  [4],  the  system  described  by  (3.1-9)  is  also  asymptotically  stable  for 


small  values  of  6jc  and  6u^,  6u? 
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3.2  The  near-optimality  characteristics. 

In  order  to  prove  that  the  suboptimal  controller  is 
really  suboptimal,  we  used  the  steepest  descent  method  to 
find  the  optimal  control  sequence  of  the  heat  pump  with  an 
initial  trial  sequence  which  is  the  same  as  that  of  the 
suboptimal  controller.  Results  of  the  search  (in  Fig. 3.1) 
show  that  the  suboptima  1  controller  sequence  is  very  close 
to  that  of  the  optimal  controller,  and  the  costs  associated 
respectively  with  these  two  controllers  are  also  close  (the 
period  involved  is  one  day). 


3.3  Simulation  and  design  methods. 

Simulation  has  become  a  powerful  tool  in  developing  as 
well  as  designing  solar  processes.  TRNSYS  program  has  been 
widely  used  in  obtaining  information  for  various  solar 
systems.  Programs  like  f-chart  have  been  used  by 
manufacturers,  and  others  to  design  solar  heating  and 
cooling  systems.  The  reason  why  simulation  is  so  popular  in 
exploring  solar  processes  is  that  physical  experiments  are 
expensive,  time-consuming,  and  generally  not  repeatable. 
However,  simulation  can  provide  us  with  much  information 
regarding  response  and  design  aspects  of  the  system.  As  John 
Duf  f i e  [ 1 ]  put  i t : 

"Properly  formulated,  simulations  can  provide  much  of 
the  same  thermal  performance  information  as  physical 
experiments  and  require  orders  of  magnitude  less  time 
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and  expense.  Effects  of  process  design  variables  can  be 
studied  systemmat ically .  Thermal  performance,  with  cost 
information,  allows  determination  of  least  cost  systems 
when  process  design  parameters  are  varied." 


However,  there  are  some  problems  confronting  the  simulator 
of  a  solar  process.  For  instance,  one  has  to  ask  what  level 
of  detail  is  appropriate  for  a  simulation?  What  kind  of 
collector  model  should  be  used?  If  more  than  necessary 
details  are  included,  computer  time  is  wasted.  On  the  other 
hand,  information  obtained  could  be  misleading  for  lack  of 
accuracy  in  modelling  the  system  components.  Meteorological 
data  is  also  a  problem,  as  Duffie  put  it: 


"...The  limited  number  of  locations  also  leaves  many 
simulators  without  data  for  locations  of  interest  to 
them,  and  there  are  still  uncertainties  in  the 
calculation  of  radiation  on  tilted  surfaces  from  data 
on  horizontal  surfaces." 


Duffie  summarized  his  remarks  concerning  simulations  as 
follows : 


"We  have  at  our  deposal  increasingly  powerful  and 
useful  tools  for  understanding  and  designing  solar 
process  systems.  Simulations  can  produce  for  us 
information  which  can  not  be  obtained  in  any  other 
pratical  way. 

The  users  of  simulation  programs,  however,  must  be 
aware  of  both  the  capabilities  and  the  limitations  of 
the  programs  they  are  using.  Misuse  of  simulation  can 
produce  more  information  on  system  in  less  time  than 
any  other  method!  Users  should  excercise  their 
engineering  judgement  in  the  use  of  simulation  based  on 
the  best  possible  knowledge  of  the  science  and 
engineering  of  the  solar  processes  they  are 
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simulating." 


The  simulation  is  done  with  the  set  of  weather  data  for 
Edmonton  from  the  first  of  October,  1967  to  the  first  of 
April,  1968.  These  data  which  are  collected  by  Environment 
Canada,  are  available  only  in  hourly  measurements; 
therefore,  in  our  simulation,  we  assume  the  weather  data  are 
unchanged  during  each  hour.  We  also  assume  that  the  average 
ambient  temperature  of  the  day  is  known  one  day  in  advance. 

A  listing  of  the  simulation  program  is  given  in  Appendix 
III.  There  are  six  subroutines  in  this  program.  To  simulate 
the  suboptimal  approach,  we  need  five  subroutines,  which  are 
READIN,  WEATHR,  A VERA ,  OPSET,  and  SYSANA.  The  WATSUN 
approach  is  simulated  through  the  use  of  subroutines  READIN, 
W7EATHR,  and  SYSWAT.  Subroutines  SYSWAT  and  SYSANA  have  their 
own  subroutines  TEMPI  and  TEMP,  respectively. 

Subroutine  READIN  is  used  to  read  in  the  parameters  of 
the  system.  Subroutine  WEATHR,  which  is  taken  from  the 
WATSUN  2  program,  is  used  to  process  meteorological  data,  as 
explained  in  the  Appendix  I.  AVERA  is  used  to  calculate  the 
daily  average  temperatures  that  are  to  be  used  in  evaluating 
the  adaptive  preview  suboptimal  controller  (see  section 
2.3.2).  Subroutine  OPTSET  computes  the  optimal  steady-state 
temperatures  of  the  high-temperature  tank,  these  values  will 
be  used  as  the  set  points  in  the  suboptimal  controller. 
SYSANA,  with  the  help  of  TEMP,  carries  out  the  simulation 
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with  information  from  subroutines  WEATHR,  A VERA ,  and  OPTSET. 
The  simulation  of  the  WATSUN  approach  is  done  by  subroutine 
SYSWAT. 

The  Euler  integration  method  is  used  with  the  sampling 
time  of  15  minutes,  or  4  samples  per  hour.  Comparing  the 
bulding  heat  loads  calculated  respectively  for  the  case  of 
4,  5,  6  samples  per  hour  (table  3.1),  a  small  error  of  0.7% 
is  reported.  This  justifies  the  use  of  4  samples  per  hour. 


3.4  Results  of  the  WATSUN  approach. 

The  strategy  of  the  WATSUN  approach  is  as  follows[2]: 

1) Check  whether  or  not  collection  is  possible.  If  no,  go  to 
the  next  step.  If  yes,  is  the  storage  tank  less  than  a 
preset  maximum  temperature?  If  yes,  tank  SI  is  provided  with 
solar  collected  energy,  if  no,  the  collected  energy  is 
dumped . 

2) Check  whether  or  not  heating  demand  exists?  If  no,  go  to 
the  next  step;  If  yes, can  tank  SI  meet  all  the  demand? 

-If  yes,  SI  supplies  heat  to  the  load  and  go  to  the 
next  step. 

-If  no,  can  tanks  SI  and  S2  meet  all  the  demand?  If  no, 
auxiliary  heat  is  used  to  supplement  SI  and  S2. 

3) Can  the  heat  pump  work  with  the  existing  temperatures  of 
tanks  SI  and  S2.  If  no,  return  to  step  1;  If  yes,  heat  pump 
takes  heat  from  tank  SI  and  supplies  it  to  tank  S2. 

The  heat  pump  used  here  is  a  3-ton  Carrier  heat  pump,  as 
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mentioned  earlier  (chapter  2).  For  lack  of  specific  data,  we 
assume  that  there  is  no  temperature  limit  on  the  heat  pump's 
range  of  operation  except  with  Tmaxas  an  upper  limt.  The 
heat  taken  from  tank  SI  to  supply  to  S2  (step  3)  is 
approximated  by  equation  (2.2-12)  with  u2,  the  heat-pump 
electrical  input,  approximated  by  the  following  equation 
( see  Fig . 3 . 2 ) : 


u  2  =  2 .  +  ( T  -T  +44 . ) /2  5 .  kW 


(3.4-1  ) 


and 


max 


3.5  Results  of  the  suboptimal  controller. 

To  use  the  suboptimal  controller  approach,  we  assume 
that  the  heat  pump  electrical  input,  u2,  can  be  changed  at 
will. 

The  strategy  for  the  control  of  the  collector  pump, 
heat  pump  and  building  circulation  pump  is  as  follows: 

1) The  collector  control  is  bang-bang,  i.e.,  whenever  there 
is  more  energy  to  be  collected  than  to  be  lost  to  the 
surroundings,  the  collector  pump  is  turned  on  to  the 
maximum,  otherwise,  it  is  switched  off. 

2) The  heat  pump  controller  is  based  upon  the  suboptimal 
controller  sequence  which  was  mentioned  earlier. 

3) The  building  circulation  pump  is  used  whenever  necessary 
to  keep  the  building  temperature  at  a  pre-determi nea  value. 
When  the  heating  load  is  not  fully  met  by  the  storage  tank, 
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auxiliary  heat  will  be  used. 

4 ) I f  the  low-side  tank  temperature  is  above  a  preset  upper 
limit,  it  indicates  a  mild  weather  condition.  Under  this 
special  condition,  the  heat  pump  will  be  turned  off  and  a 
mass  exchange  between  the  two  tanks  instituted,  thus 
effectively  making  the  system  a  one-tank  pure  solar  system. 
The  simulation  period  is  150  days  long.  It  is  from  October 
1,  1967  to  April  1,  1968.  This  period  can  be  roughly  divided 

into  three  smaller  periods.  The  first  period  is  from  the  1st 
day  to  the  80th  day.  The  second  period  is  from  the  81th  day 
to  the  130th  day,  and  the  third  period  includes  days  of  the 
rest  of  the  simulation  period  (see  Fig. 3. 3  and  Fig. 3. 4).  The 
simulation  results  of  the  WATSUN  and  the  suboptimal 
approaches  are  plotted  in  Fig. 3. 5  and  3.6  which  show  that 
the  energy  consumed  by  the  WATSUN  approach  is  more  than  that 
of  the  suboptimal  controller  for  the  first  period,  and  for 
the  third  period.  In  the  first  period,  the  weather  is 
considered  good  by  the  suboptimal  controller,  hence, 
strategy  4  is  taken.  Therefore,  in  this  period,  the  saving 
of  purchased  energy  is  mainly  due  to  the  doubling  of  the 
mass  storage  of  the  suboptimal  controller.  In  the  following, 
we  will  look  at  the  performances  of  the  two  approaches  only 
in  the  second  and  the  third  period  (Fig. 3. 7  and  Fig. 3. 8)  to 
see  whether  the  suboptimal  controller  consumes  more  or  less 
energy  in  less  favourable  conditions  of  weather.  In  the 
third  period  of  interest,  130th  day  to  150th  day,  the 
weather  is  relatively  better  than  the  second  period  (but 
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worse  than  the  first),  and  the  suboptimal  controller  (using 
strategy  2)  makes  a  saving  of  10%  in  purchased  energy  ,as 
seen  in  Fig.3.7(the  total  saving  for  the  entire  time,  150 
days,  is  31%).  Fig. 3. 7  also  shows  that  in  the  cold  period, 
80th  day  to  130th  day,  the  suboptimal  controller  and  the 
WATSUN  approach  consume  almost  the  same  amount  of  purchased 
energy.  In  other  words,  there  is  no  benefit  using  the 
suboptimal  controller  when  the  weather  is  so  bad.  If  we 
increase  the  collector  area  from  100  to  300  m2,  using 
strategy  2,  37%  saving  of  purchased  energy  is  obtained 
(Fig. 3. 8).  But  most  of  these  savings  come  from  the  third 
period.  This  ,once  more,  confirms  the  above  conclusion  that 
unless  the  weather  is  very  bad  the  suboptimal  controller 
always  outperforms  the  WATSUN  algorithm.  It  should  also  be 
noted  that  the  second  period  is  for  the  months  of  December 
and  January  which,  for  Edmonton,  are  not  likely  to  be 
suitable  for  solar  heating  in  any  case. 

Table  3 . 1 
samples  per  hour 
4  5  6 


Building  load(GJ)  69.39  69.34 


69.33 


CONTROL  VALUE 
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TIME  (HOUR) 


Fig.3.1:  Suboptimal  and  Optimal  Controller  Values 
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Fig. 3. 2:  Heat  Pump  Electrical  Input 


64 


Fig. 3. 3:  Average  Insolation 
(per  ten  days) 


(per  ten  days) 


PURCHASED  ENERGY  ( MJ  )  X 1 O4 


65 


Fig. 3. 5:  Simulation  results 
(collector  area=100  m2) 


PURCHASED  ENERGY  C  MJ  )  X 1 04 


66 


Fig. 3.6:  Simulation  results 
(cummulative  energy  cost,  collector  area=100  m2  ) 


PURCHASED  ENERGY  ( MJ )  X 1 04 
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Fig. 3. 7:  Simulation  results 
(2nd  and  3rd  periods,  collector  area=100  m2) 


PURCHASED  ENERGY  ( MJ  ) 
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Fig. 3. 8:  Simulation  Results 
(collector  area=300  m2 ) 
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4.  THERMAL  and  ECONOMIC  ANALYSIS 


4.1  Thermal  sensitivity  analysis. 

The  performance  of  a  solar  heating  system  depends  on 
many  factors,  some  of  which  are  beyond  our  control,  e.g.  the 
ambient  temperature  or  the  insolation.  But  there  are 
possibilities  to  improve  the  performance  of  a  system.  For 
example,  it  has  been  shown  that  the  system  efficiency  can  be 
enhanced  with  the  proper  sizing  of  system  components.  The 
purpose  of  this  study  is  to  identify  areas  for  possible 
performance  improvement  and  (or)  cost  reduction  through 
utilization  of  thermal  and  economic  models. 

The  saving  percentage  is  defined  as  the  percentage  of  the 
system  heat  load  which  is  provided  by  the  solar  energy;  or 
spec i f i ca 1 ly ,  it  is  equal  to  , in  percent,  the  difference 
between  the  total  building  heat  load  and  the  purchased 
energy,  divided  by  the  total  building  heat  load. 


4.1.1  Mass  storage . 

For  a  given  collector  area,  various  values  of  the 
storage  mass  per  unit  collector  area  were  used  to  calculate 
the  savings  percentage  (it  is  assumed  that  the  tanks  are 
always  of  equal  mass).  Using  simulations  as  described  in  the 
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previous  chapter,  with  the  same  particular  set  of  weather 
data,  each  value  of  the  storage  mass  was  used  to  find  the 
corresponding  savings  percentage  (Fig.  4.1).  In  Fig. 4.1,  we 
note  that  when  the  collector  area  is  small,  the  system 
performance  becomes  worse  as  the  storage  mass  per  unit 
collector  area  increases  from  40.  to  100  Kg  per  one  square 
meter  of  collector  area.  The  phenomenon  is  reversed  when  the 
collector  area  is  large.  This  is  due  to  the  fact  that  if  the 
tank  storage  mass  is  larger,  the  average  tank  temperature 
becomes  lower,  and  this  leads  to  a  worse  heat-pump  COP,  but 
to  a  better  collector  efficiency.  Hence,  it  follows  that 
better  performance  of  the  collector  does  not  compensate  for 
a  worse  COP  in  the  case  of  a  small  collector.  The  reverse  is 
true  for  a. large  collector,  because  if  a  collector  is  large 
enough,  solar  collection  would  out-perform  the  heat  pump[ 1 ] . 


4.1.2  Collector  area. 

For  a  given  storage  capacity  per  unit  collector  area, 
the  savings  percentage  is  directly  proport ional  to  the 
collector  area  (Fig4.2).  A  larger  collector  helps  capture 
more  solar  energy,  and  therefore  more  energy  can  be  stored 


to  be  sent  to  the  load  later. 
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4.1.3  Heat  exchanger  coefficient. 

The  heat  exchanger  coefficient,  which  is  positive  and 
always  less  than  or  equal  to  unity,  can  cause  a  difference 
of  25%  in  the  saving  percentage  if  the  coefficient  is 
increased  to  1.  from  its  lower  value  of  0.1  (Fig. 4. 3).  The 
heat  exchanger  is  used  to  transfer  heat  from  the  liquid 
inside  the  duct  to  the  room  airspace.  A  higher  heat 
exchanger  coefficient  means  a  larger  amount  of  heat  to  be 
t ransf errable  and  a  less  amount  of  auxiliary  heat  to  be 
needed . 


4.1.4  Collector  covers. 

Using  the  numbers  in  [5],  the  savings  results  of  the 
system  with  zero-,  single-,  and  doubled-glazed  are  tabulated 
in  table  4.1.  The  zero-glazed  is  always  out -pe r f ormed  by  the 
single-,  and  the  doubled-glazed.  For  large  collector,  the 
double-  and  the  single-glazed  save  almost  the  same  amount  of 
energy,  but  for  smaller  collector,  the  single-glazed  can 
cause  a  difference  of  0.16%  in  saving  of  purchased  energy, 
therefore,  the  single-glazed  colector  is  the  most 
appropriate.  This  agrees  with  the  conclusion  of  [6]. 


4.1.5  Building  loop  pump  capacity. 
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Supposing  that  we  encounter  a  scenario  as  follows: 
there  is  a  large  instantaneous  heating  demand,  which  the 
high-side  tank  is  capable  of  meeting,  but  the  building  loop 
pump  capacity  is  not  sufficient  to  transfer  fully  to  the 
room  the  amount  of  energy  needed.  In  that  case,  auxiliary 
heat  must  be  used  to  keep  the  room  at  desired  temperature, 
thus  lowering  the  saving  percentage,  as  can  be  seen  in 
Fig. 4. 4.  However,  the  building  loop  pump  only  aids  in 
transferring  heat  to  the  load  on  the  condition  that  the  tank 
stores  enough  energy  and  a  heating  demand  exists;  therefore, 
there  is  a  limit  on  the  increase  of  the  saving  percentage 
when  the  loop  pump  capacity  is  increased.  This  explains  the 
saturation  branch  of  the  curve  in  Fig. 4. 4.  Approximately  16% 
more  of  saving  percentage  is  recorded  here.  An  effective 
collector  area  of  20  m2 ,  a  mass  storage  of  75  Kg  per  unit 
collector  area  (m2),  and  a  heat  exchange  coefficient  of  1. 
were  used  to  plot  fig. 4. 4. 


4.1.6  Tanks  heat  loss  coefficient. 

The  change  of  the  tanks  heat  loss  coef f ic ient ( a 
and  a^  )  does  not  cause  a  significant  change  in  the 
performance  of  the  system.  Less  than  2%  difference  of  the 
saving  energy  is  found  when  the  heat  loss  coefficients  are 
changed  from  5  to  15KJ/hr°C  (Fig. 4. 5). 
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4.1.7  Tilt  slope . 

A  variation  of  the  tilt  slope  could  result  in  a  worse 
or  a  better  saving  percentage.  In  Fig. 4. 6,  the  optimum  slope 
is  found  to  be  in  the  neighborhood  of  the  region  lattitude. 
This  agrees  with  the  rule  of  thumb  that  the  optimum  tilt 
slope  is  in  the  range  from  ( lat i tude- 1 0 °  )  to 
( lat i tude+ 10°  ) . 


4.1.8  Heat  pump  capacity. 

The  optimum  value  of  the  heat  pump  capacity  is  around 
5kW  (Fig. 4. 7).  The  savings  percentage  is  slightly  lower  if 
the  heat  pump  capacity  is  increased  beyond  5  KW.  This  is 
understandable  if  we  we  note  that  in  our  simulation  time, 
there  are  fourteen  days,  the  93rd  day  to  the  106th  day, 
without  sunlight.  Hence,  in  this  period,  the  steady-state 
optimal  temperatures  of  the  two  tanks  do  not  exist  (Chapter 
II),  though  we  can  still  calculate  the  pseudo-optimal  values 
of  the  high-temperature  tank  by  equation  (2.3-45). 

Therefore,  any  successful  efforts  by  the  heat  pump  to  reach 
these  pseudo-optimal  values  will  result  in  more  consumption 
of  purchased  energy.  When  the  heat  pump  capacity  is  less 
than  5  KW,  the  system  will  not  be  able  to  reach  the 
pseudo-optimal  points.  The  reverse  is  true  for  systems  with 
higher  heat  pump  capacity.  However,  if  we  shut  off  the  heat 
pump  completely  during  the  cold  period  this  will  degrade  the 


0 
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system  performance.  The  reason  is  that  the  heat  pump  will 
work  with  bad  values  of  the  COP  (because  of  low  values  of 
the  low-side  temperature)  when  the  weather  becomes  milder, 
thus  decreasing  the  saving  percentage.  This  is  another  area 
where  further  work  is  needed  to  improve  the  performance  of 
the  suboptimal  controller. 


4.2  Economic  analysis. 

The  thermal  performance  of  the  system  was  dealt  with  in 
the  preceding  discussion,  but  the  real  impact  on  the  part  of 
the  designer  lies  on  the  question  of  cost.  No  matter  how 
good  the  system  performance  is,  it  will  not  be  realizable  if 
the  cost  is  unpractically  high.  According  to  J.Duffie  and 
W . A . Bee kman [ 2 ] ,  the  annual  cost  of  a  solar  heating  system 
consists  of  many  elements  such  as  the  annual  cost  of 
ownership,  the  annual  cost  of  operation,  and  the  yearly  cost 
of  maintenance.  As  Duffie  and  Beckman  put  it: 


"The  major  annual  cost  of  a  solar  heating  system, 
without  auxiliary  energy,  includes:  the  annual  cost  of 
operating  the  system;  the  cost  of  power  for  the  pumps, 
blowers,  and  so  on;  and  the  yearly  cost  of  maintenance. 
The  annual  cost  of  ownership  includes  cost  associated 
with  the  initial  investment,  that  is,  interest  on  the 
investment  and  its  repayment  over  a  specified  number  of 
years  related  to  its  lifetime.  The  sum  of  these  is 
usually  taken  as  a  fixed  percentage  of  investment  each 
year;  for  example,  for  a  twenty  years  amortization  and 
8%  interest  rate,  the  annual  cost  is  0.10185  of  the 
investment . 

Operating  costs  are  primarily  for  power  requirement  for 
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pumping  water  and  moving  air  in  the  system,  summed  over 
the  yearly  operating  time  of  the  system.  Maintenance 
costs  include  repairs,  replacement  of  glass  in 
collectors,  or  any  other  costs  of  keeping  the  system  in 
operating  condition.  Consideration  of  these  costs  leads 
to  the  conclusion  that  maintenance  must  be  minimized  if 
solar  heating  is  to  be  economically  viable, 
particularly  when  labour  costs  are  to  be  charged  as 
part  of  the  maintenance  expense." 


Similarly,  the  cost  of  auxiliary  also  consists  of  the 
equipment  cost  (e.g.  furnace),  the  fuel  cost,etc... 

There  are  other  factors  that  can  further  complicate  the 
matter,  for  example,  insurance  or  real  estate  taxes  or  tax 
relief  on  residential  use  of  solar  energy.  For  lack  of  data, 
and  for  simplicity,  the  cost  analysis  of  a  solar  heating 
system,  which  is  based  on  [3],  and  [4]  can  be  broken  down  as 
follows  (Canadian  dollar  is  used): 

-The  cost  of  the  heat  pump  is  $3000.  The  annual  operating 
cost  cost  of  the  heat  pump  is: 

C,=QhpxCel 

where , 

Q  =Amount  of  electricity  consumed  by  heat  pump 

hp 

[KWhr ] 

Cei=Utility  energy  rate  [$/KWhr] 

-The  collector  cost: 

C2=A  xC  , 

c  cl 

where , 

A  =  Effective  collector  area  [m2 ] 

Cci  =  Collector  cost  [$/m2] 


-The  storage  cost: 
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=  Storage  cost  [$/Kg] 

Ms=  Storage  mass  [Kg] 

-The  cost  of  pumps,  valves,  pipes,  heat  exchanger 

C4=200.+10xAc 

-The  controller  cost: 


C  5  = 1 50 . 

-The  capital  cost  of  the  electric  heater  (used  as 
auxiliary  heater)  is  $200. 

-The  cost  of  auxiliary  energy  used: 


C‘  QauxxCel 

If  I  is  the  capital  cost  of  the  system,  for  given 
period  of  N  years,  and  a  given  annual  interest  of 

amortized  capital  cost  in  year  j  i  s [ 4 ] : 

N  N 

il  (  1  +  i)  /[ ( 1 .+i)  -1 . ] 


amort i ze  d 
i % ,  the 


For  an  annual  real  inflation  rate  of  fuel  at  r%,  the 
operating  cost  in  year  j  i  s [ 4  ]  : 

J 

(Qa  +Q,  )Cpl  (  1  .+r) 
aux  hp  el 

So,  the  total  cost  of  the  system  in  year  j  is: 

C  .  = (  3000 . +  A  xC  n  +  C_,  M  +200  .  +  1  OxA  + 

SI  C  Cl  St  s  c 

N  N  J 

1  50  . +200 . )xi ( 1  +  i )/[ ( 1  +  i )- 1 ]  +  (Q  +QU  )Cpl(1+r) 

aux  hp  ei 


(4.2- 


The  total  cummulative  cost  of  the  solar  system  is: 

N 

C  =tc  . 
s  4^  sj 

1=1 

Assuming  that  the  mortgage  life  ,N,  is  equal  to  the  system 
age,  the  annual  interest  rate  is  12%,  and  the  annual  real 
inflation  rate  of  electricity  and  natural  gas  is  10%.  Using 
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the  values  in  [3], we  have: 

Cc-|  =$60/m2  , 

C$t=$0. 132/Kg, 

Cg-j  =$0 . 00335/KWhr  .  (Edmonton  Power  price) 

The  assumed  conventional  system  uses  natural  gas  (which  is 
widely  available  in  Alberta)  to  heat  the  house.  Suppose  the 
cost  of  the  furnace  is  $1000,  which  is  also  financed  in  the 
same  amortized  period  and  at  the  same  interest  rate  as  that 
of  the  solar  system. 

Similarly,  the  total  cost  of  the  conventional  system  at  year 
j  is: 

N  N  j 

C  .  =  1 OOO.xi ( 1  +  i )/[ ( 1  +  i )  -  1 ] +0  Cn  (1+r) 

cvj  ~cv  9 

The  price  of  natural  gas  in  Edmonton,  Cg  ,  is  $2.94/GJ 
(Northwestern  Utilities  price,  tax  included). 


4.2.1  The  viability  of  the  solar  system. 

A  computer  program  is  written  to  compare  the  costs 
associated  respectively  with  the  above  two  systems  when  the 
solar  system  age  varies  from  10  to  30  years.  With  the 
collector  area  of  20  m2,  and  the  storage  mass  per  unit 
collector  area  of  75  Kg,  the  results  are  calculated  and 
plotted  in  Fig. 4. 8,  it  is  noted  that  from  this  figure  that 
the  solar  heating  system  is  not  economically  viable.  If  the 
collector  area  is  increased  to  40  m2 ,  the  solar  heating 
system  is  not  viable  either,  as  seen  in  Fig. 4. 9.  Hence,  the 


’ 
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economic  viability  of  solar  heating  is  possible  only  if  the 
solar  industry  is  developed  enough  to  decrease  the 
investment  cost  on  solar  equipments. 


4.2.2  Costs  versus  mass  storage  and  collector. 

Using  the  same  amortized  period,  annual  interest  and 
real  inflation  rate,  the  system  cummulative  costs  change 
with  respect  to  respective  variation  of  the  collector  area 
and  the  mass  storage  is  plotted  in  Fig, 4. 10  and  Fig. 4. 11. 
These  figures  show  that  the  collector  and  the  mass  storage 
are  dominating  factors  in  the  total  system  costs.  Therefore, 
we  will  not  be  able  to  base  our  design  criterion s  on  the 
energy-saving  aspect  alone,  unless  at  some  future  time  the 
capital  costs  spent  on  the  collector  and  the  mass  storage 
are  down. 


4.2.3  Costs  versus  Heat  Pump  Capacity 

Based  upon  data  provided  by  AAF  Ltd.  (5718 
Edmonton,  Alberta),  the  heat  pump  price  per  kW 
approximated  by  (Fig. 4. 12): 

P  =250 ( W  ) +900 

hp  max 

where  P,  =price  of  heat  pump  [$C], 

hp 


103  St  . 
i  s 


(4. 2. 3-1 ) 
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W  =maximum  electrical  input  to  heat  pump 
max  ^  t'  t' 

[kw] 

Replacing  3000.  in  equation  (4.2-1)  by  P  ,  and  with  the 

hp 

same  assumption  regarding  interest  rate,  amortized  period, 
etc.,  the  cummulative  system  costs  are  calculated  with 
various  values  of  W  x  and  plotted  in  Fig. 4. 14  (collector 
area  =  20  m2  )  and  Fig. 4. 13  (collector  area=100  m2).  Like 
collector  and  mass  storage,  the  heat  pump  cost  also  plays  a 
dominant  role  in  the  system  cost  . 


4.3  Recommendation. 

The  preceding  discussion  dealt  with  the  thermal  and 
economic  aspects  of  a  solar  heating  system.  Based  on  these 
results,  a  few  conclusions  can  be  drawn.  However,  it  should 
be  noted  here  that  due  to  the  unavailability  of  suitable 
data(  for  example  the  only  complete  set  of  weather  data 
available  to  us  is  from  October  1,  1967  to  April  1,  1968), 

the  previous  results  as  well  as  the  following  comments 
should  be  taken  with  caution. 

-The  tilt  slope  of  the  collector  should  be  made  as  close  as 
possible  to  the  lattitude. 

-Single-glazed  collector  should  be  used. 

-The  building  loop  pump  capacity  should  be  increased  to  the 
neighborhood  of  a  proper  value  if  this  does  not  lead  to  an 
increase  in  total  cost.  Here,  for  lack  of  specific  pump 
cost,  it  is  not  possible  to  investigate  the  economic 
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-The  heat  exchanger  coefficient  should  be  kept  as  close  to 
1.  as  practically  and  economically  possible.  The  economic 
aspect  of  the  system  response  due  to  variation  of  the  heat 
exchanger  is  not  studied  for  lack  of  data. 

-The  heat  loss  coefficient  of  the  tanks  can  cause  a  2% 
decline  in  the  saving  percentage  for  a  change  from  5  to 
12KJ/hr  C.  This  points  out  the  need  for  a  better  insulated 
tank  though  the  saving  is  not  much. 

-The  increased  storage  mass,  heat  pump  capacity,  and 
collector  area  respectively  help  boost  the  saving  percentage 
but  also  make  the  system  more  costly. 

-The  use  of  solar  energy  in  Edmonton  is  not  viable. 

In  fact,  there  are  many  other  ignored  factors  that  can 
further  complicate  the  matter  and  enhance  or  destroy  the 
economic  aspect  of  the  system.  They  are  such  as  tax 
incentive,  maintenance  costs,  property  tax  ,  e tc  .  . . and  they 
are  beyond  the  scope  of  this  thesis. 


4.4  Further  research. 

As  mentioned  earlier  (section  2.4.2  and  4.1.8),  the 
suboptimal  controller  is  not  "smart"  enough  to  take  in  to 
account  the  existence  condition  of  the  steady-state  optimal 
tank  temperatures.  However,  we  believe  that  there  is  the 
possibility  to  upgrade  the  thermal  performance  of  the  system 
if  further  work  is  done  on  this  area. 
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Thus  far  we  assume  that  the  process  involved  is 
deterministic  in  nature.  In  reality,  the  weather  forecast 
data,  which  are  only  probable  values  are  subject  to  random 
effects  and  it  is  essential  that  these  effects  are  taken 
into  account  in  our  calculations.  The  stochastic  optimum 
theory  may  be  useful  in  studying  the  stochastic  nature  of 
the  process. 

Another  problem  which  remains  unsolved  is  the 
implementation  of  the  suboptimal  controller  on  a 
microprocessor.  The  calculation  of  the  controller  values 
(equation  (2. 4. 2-1))  is  relatively  simple.  However,  if  we 
want  to  avoid  developing  arithmetic  subroutines  (i.e. 
addition,  subtraction,  multiplication,  etc.),  a 
number-oriented  microprocessor  (calculator  chip)  can  be  used 
as  in  [7],  According  to  T.B.Kent  et  al[8],  as  far  as 
technical  problems  to  design  a  microprocessor  controller  are 
concerned  they  are  non-existent.  In  other  words,  the 
implementation  of  this  suboptimal  controller  is  certainly 
possible.  The  only  problems  left  are  to  develop  proper 
hardware  and  software  to  support  it. 

Since  utility  companies  might  offer  a  lower  price  on 
electricity  produced  during  off-peak  hours,  it  would  be 
economically  advantageous  to  use  the  auxiliary  electrical 
energy  for  the  system  during  that  period  while  reducing  peak 
hour  demands  of  electricity.  The  study  of  this  problem  can 
also  be  extended  to  examine  the  impact  of  the  system 
electrical  load  on  the  power  network. 
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The  above  are  but  a  few  related  areas  that  needs  more 
studies . 


4 . 5  Conclus i on . 

Using  the  KUHN-TUCKER  theorem,  we  have  obtained  the 
expression  of  the  steady-state  optimal  values  of  the 
high-side  temperature,  the  low-side  temperature,  and  the 
heat-pump  controller.  As  weather  conditions  change,  the 
optimal  heat  pump  controller  sequence  is  numerically 
evaluated.  Based  upon  this  result  and  upon  the  results  of  a 
process  of  trial  and  error,  a  suboptimal  adaptive  preview 
controller  is  proposed.  Its  simulation  result  is  analysed 
and  compared  with  that  of  the  WATSUN  approach.  Savings  of 
31%  or  more  are  recorded  for  the  suboptimal  controller.  The 
thermal  and  economic  behaviour  of  the  system  are  also 
investigated  to  identify  areas  for  possible  improvement  of 
the  system  performance.  Further  research  is  also  proposed. 
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Table  4 . 1 


saving 

percentage 

Collector  area 

100  m2 

20  m2 

Zero-glazed 

36.61% 

30.80% 

S ingle-glazed 

42.83% 

30.86% 

Double-glazed 

42.84% 

30 . 72% 

SAVING-  PERCENTAGE 
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Fig  4.1:  Saving  percentage  versus  mass  storage 
(collector  area= 1 00 , 50 , 20  m2  ) 
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Fig  4.2:  Saving  percentage  versus  collector  area 
(collector  area=  20-120  m2,  mass  storage=  75  kG/m2) 
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Fig  4.3:  Saving  percentage  versus  heat  exchanger  coefficient 

(collector  area=20  m2 ) 
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Fig  4.4:  Saving  percentage  versus  building  loop  pump 

capacity 

(collector  area=20  m2 ) 


SAVING-  PERCENTAGE 
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HEAT  LOSS  COEE.  (AC) 


Fig  4.5:  Saving  percentage  versus  tanks  heat  loss 

coef  f ic ient (a  ,  ,  a  ) 

h  c 

(collector  area=20  m2 ) 


SAVING  PERCENTAGE 
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Fig  4.6:  Saving  percentage  versus  tilt  slope 
(collector  area=20  m2 ) 


SAVING  PERCENTAGE 
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HEAT  PUMP  CAPACITY  (  KW  ) 


Fig  4.7:  Saving  percentage  versus  heat  pump  capacity 

(collector  area=  100  m2 ) 


s  0  l  X  ($3)  A3N0U 
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Fig. 4. 8:  Cummulative  costs  of  the  solar  and  conventional 


systems . 

(collector  area=20  m2) 


MONEY  (C$)  X105 
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Fig. 4. 9:  Cummulat i ve  costs  of  the  solar  and  the  conventional 


systems . 

(collector  area=40  m2) 


CUMMULRT I VE  COST  ( C$ )  X104 
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Fig  4.10:  Cost  versus  mass  s tcrage ( Kg/m2 ) . 
(collector  area= 1 00 , 50 , 20  m2 ) 


CUMMULRT I VE  COST  ( C$ )  X 1 04 
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Fig  4.11:Cost  versus  collector  area, 
(mass  storage= 1 00 , 75 , 50  kg/m2) 
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Fig. 4. 12:  Heat  pump  price  versus  capacity 
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H.p.  CnPRCITY  (KN) 


Fig. 4.  13:  cummiative  cost  versus  heat  pump  capacities 

(collector  area=100  m2 ) 
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Fig. 4.  14:  Cummulative  cost  versus  heat  pump  capacities 

(collector  area=20  m2) 
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APPENDIX  I 


Meteorological  data  processing 


The  hourly  measurements  of  the  radiation  and  the  ambient 
temperature  are  needed  for  using  the  suboptimal  controller. 
It  is  the  radiation  incident  on  the  tilted  collector  that  is 
needed  for  calculating  the  collected  energy  (and  the  gain  of 
the  controller),  but  the  only  values  provided  by  Environment 
Canada  are  the  total  radiation  (H)  incident  on  a  horizontal 
surface.  A  program  based  upon  the  following  algorithm  was 
written  to  process  these  data[1]. 

First,  calculate  the  beam  (H  )  and  the  diffused(H  ) 

b  d 

components  of  H: 


and , 


where,  K  can  be  calculated  by  using  a  correlation: 


K=  0 . 177 


if 


Ky>0 .75 
0.35<Kt<0.75 
C  .  <Kj<0 . 35 


K=  1  .  56-  1  .  84K.J.  i  f 
K=  1  . -0 . 249K  T  if 


and  Kt  =H/H0 

with  H Q  being  the  instantaneous 


rad i a  t i on 


on  a  horizontal  surface: 


H  Q=Sc [  1+0.033cos(2ttN/365)  ]  x  [  coscocos</?cos<5  +  sinps  in  <5  ] 
where,  Sc=  solar  constant 
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1  10 


N=  day  of  the  year 
<5  =  solar  declination 
v  =  lat t i tude 
co  =  hour  angle 

The  H.J.  can  be  determined  from: 

H  =R  Hl  +( 1+cosS)xH  /2.+  p ( 1-cosS)H/2. 

1  b  b  d 

S=  tilted  angle  of  collector 
p=  ground  reflectance  factor 

R  is  the  ratio  of  beam  radiation  on  a  tilted  surface  to 

b 

that  on  a  horizontal  surface;  and  if  the  collector  is  facing 
towards  the  equator,  R  is  calculated  as  follows: 

R  =[cos(f  -  s)cos5cosa>  +sin(<P-s)sin  h  3/ 

b 

[ cos^cos 5cos  co  +sin<psinS] 
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APPENDIX  II 


Program  to  find  the  optimal  controller  sequence 


The  following  program  is  originally  form  [1],  with  weather 
data  and  system  parameters  being  purely  hypothetical  [2]. 

C  PROGRAM  TO  OPTIMIZE  CONTROL  FOR 

C  HEAT  PUMP  SYSTEM 

DIMENSION  U1 (51 ) , U2(51 ) , TH ( 5 1 ) , TL ( 5 1 ) , 

1  RL ( 5  1  )  , COST (51), PCOST (51), GRA (51), 

2  U2L (51), CO STL (51 ) ,DIR(51 ) ,DI  (51 ) , 

3  QAUX (51), GRAD (51) , RH ( 5 1 ) 

C . 

C  PART  ONE 

C  ESTABLISHMENT  OF  THE  INITIAL  TRAJECTORY 

C 

c . . 

DT=  0 .  1 5 
TB  =  2  0 . 

TA  =  - 2  0 . 

TH (  1  ) =4  0 . 

TL  ( 1 )  =  5  . 

S=  1  .  0 
I TER= 1 
J  =  0 


COST ( 1 ) =0 . 


TCMI N= 1000. 


1  12 


DO  100  1=1,50 
I F ( I .GT.  1 5 ) S  =  0 . 

I F ( I .GT.  1 5 ) TA  =  - 2  0 . 

I F ( I . GT . 35 ) S= 1 .0 
I F ( I . GT . 35 ) TA=~20 . 

DELT=TB-TA 

U1  (I )  =  0 . 5*DELT/(TH( I )-TB) 

I F (U 1 ( I )  .GT.  1 )U1  (I  )  =  1  . 

QAUX ( I ) =0 . 1*DELT-0.2*(TH(I )-TB)*U1 (I  ) 
UNUM=0.2*U1 ( I ) * ( TH ( I ) -TB ) +0 . 02* ( TH ( I ) -TA ) 
UDEN=20* ( 1 . - . 0 1 6  * ( TH ( I ) -TL ( I ) ) ) 

U2 ( I ) =UNUM/UDEN 
I F ( U  2 ( I ) .GT.  1 )U2(I )  =  1 . 

I F ( ( TH ( I )-TL(I ) ) .GT.45)U2(I )=0. 

U2L ( I ) =U2 ( I ) 

RH ( I )=UDEN*U2(I ) -UNUM 

TERM  1 =  2  0 . *U2 ( I ) * ( . 7 1 5- . 0 1 6* ( TH ( I )-TL(I ) ) ) 
+  +0.02*(TL(I )-TA) 

TERM 2= 10.*(S~0.025*(TL(l )-TA) ) 

IF (TERM2 .LT. 0 )TERM2=0 . 

RL( I ) =TERM2-TERM1 

TH  ( I  +  1  ) =TH ( I ) +DT*RH ( I ) 

TL (1+1 ) =TL ( I ) +DT*RL ( I ) 

COST (1+1 ) =COST ( I ) +  QAUX ( I )+5.7*U2(I ) 

100  CONTINUE 


TCOST=COST ( 5 1 ) 


WRI TE ( 6 , 1 ) I  TER , J , TCOST 
WRITE (6 , 4 ) 

WRI TE (6,2) (U2(I ) ,1=1 ,50) 

WRI TE ( 6 , 6 ) 

WRITE (6, 3) (TH(I ) ,1=1 ,50) 

WRI TE ( 6 , 8 ) 

WRITE (6, 3) (TL (I), 1=1, 50)  . 

C . 

C  PART  TWO 

C  EVALUATION  OF  THE  GRADIENTS 

C  AND  DOING  ONE-DI MENS I ONED  SEARCH 

C . 

GMA2= 1 . 

GM2= 1 . 

125  JP=  3  5 

GA 1 =GA2 
DMAG2=0 . 

GMAG2  =  0 . 

150  PCOST ( JP ) =COST ( JP ) 

U2 (JP)=U2L(JP)+0 .01 
DO  200  I = JP , 50 
I F ( I .GT. 1 5 ) S  =  0 . 

IF ( I  .GT.  1 5 ) TA  =  - 2  0 . 

I F ( I . GT . 3  5 ) S=  1  .0 
IF ( I .GT. 35 )TA=-20 . 

DELT=TB-TA 

U1 (I )=0 . 5*DELT/(TH( I )-TB) 
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I F ( U 1 ( I ) . GT . 1 ) U 1 ( I  )  =  1  . 

QAUX ( I ) =0 . 1*DELT-0.2*(TH(I )-TB)*U1 (I ) 
UNUM=0.2*U1 (I )*(TH(I ) -TB ) +0 . 02* ( TH ( I )-TA) 
UDEN=20* ( 1 . 0 1 6* (TH( I ) -TL ( I ) ) ) 

I F ( I . LE .35) GO  TO  175 

U2 ( I )  =  ( 1 .  +  1 . * ( 4  0 . -TH ( I ) ) ) *UNUM/UDEN 
175  CONTINUE 

I F ( I . EQ . JP ) GO  TO  180 
I F ( U2 ( I ) .GT. 1 )U2(I )=1 . 

180  IF ( (TH( I ) -TL( I ) ) .GT. 45 )U2 ( I ) =0 . 

I F ( U  2 ( I )  . LT . 0 ) U2 ( I ) =0 . 

RH ( I ) =UDEN*U2 ( I ) -UNUM 

TERM  1  =  20. *U2 ( I )  *  ( .715~.016*(TH(I )-TL(I ) ) ) 
+  +  0 . 0  2  * ( TL ( I ) -TA ) 

TERM2= 1 0 . * (S-0 . 025* ( TL ( I )-TA) ) 

IF (TERM2 .LT. 0 )TERM2=0 . 

RL ( I ) =TERM2-TERM1 
TH (1  +  1  ) =TH ( I ) +DT*RH ( I ) 

TL (1+1 ) =TL ( I ) +DT*RL ( I ) 

PCOST (1  +  1 ) =  PCOST ( I ) +  QAUX ( I )+5.7*U2(I ) 

200  CONTINUE 

GRA ( JP) =GRAD( JP) 

GRAD ( JP ) = ( PCOST ( 5 1 ) -TCOST ) /0 . 0 1 
GMAG2=GMAG2+ ( GRAD ( JP ) ) **2 
DIR( JP)=-0 . 5*GRAD( JP)+0. 5*DI ( JP) 
DMAG2=DMAG2+ ( DI R ( JP ) )**2 


DI ( JP)=DIR( JP) 
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I F ( JP.EQ. 1 ) GO  TO  250 
JP=JP- 1 
GO  TO  150 
250  CONTINUE 
GM2=GMA2 
GMA2=GMAG2 
R=GMA2/GM2 
JMI  N  = 1 

DO  600  J=  1  ,  1 0 
G=0 . 0  1  *J 

270  DO  500  1=1,50 

U2 ( I ) =U2L ( I ) +G*DI R ( I ) / ( DMAG2 ) *  *  0 . 5 
IF (I .GT. 15) S=0 . 

I F ( I .GT. 1 5 ) TA  =  - 2  0 . 

I F ( I . GT . 3  5 ) S= 1  .0 
IF ( I .GT. 35)TA=-20 . 

DELT=TB-TA 

U 1  (I ) =  0 . 5*DELT/(TH( I )-TB) 

IF(U1  (I )  .GT.  1 )U1  (I  )  =  1  . 

QAUX(I ) =0 . UDELT-0 . 2* (TH( I )~TB)*U1 (I ) 
UNUM=  0 . 2  *U 1 (I )*(TH(I ) -TB ) + 0 . 0 2 * ( TH ( I )-TA) 
UDEN  =  2  0  * (  1 . - . 0 1 6* ( TH ( I )-TL(I )  )  ) 

IF (I .LE. 35)GO  TO  275 
U2 ( I )  =  ( 1 .  +  1 . * ( 4  0 . -TH ( I ) ) ) *UNUM/UDEN 
275  CONTINUE 

IF ( U  2 (I ) .GT.  1 )U2(I )  =  1 . 

I F ( ( TH ( I ) -TL ( I ) ) .GT.45)U2(I )=0. 
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I F ( U2 ( I ) . LT . 0 ) U2 ( I  ) =0 . 

RH ( I ) =UDEN*U2 ( I ) -UNUM 

TERM1 =20 . *U2 ( I ) * ( . 7 1 5- . 0 1 6* ( TH ( I ) -TL ( I ) ) ) 
+  +0.02*(TL(I )-TA) 

TERM2=10.*(S-0.025*(TL(I  )-TA) ) 

IF (TERM2 .LT. 0 )TERM2=0 . 

RL ( I ) =TERM2-TERM1 
TH(I  +  1 ) =TH ( I ) +DT*RH ( I  ) 

TL (1  +  1 ) =TL ( I ) +DT*RL ( I  ) 

COST (1+1 )=COST(I ) +QAUX ( I )+5.7*U2(I ) 

500  CONTINUE 

TCOST=COST ( 5 1 ) 

IF (KEEP. EQ. 1 )GO  TO  610 
IF (TCOST.LT. TCMIN) JMIN=J 
IF (TCOST.LT. TCMIN)TCMIN=TCOST 
600  CONTINUE 

G  =  0 . 0 1  * JMI N 
KEEP= 1 
GO  TO  270 
610  I TER= I TER+ 1 

DO  620  K= 1  ,  30 
620  U2L ( K ) =U2 ( K ) 

KEY=ITER/50 
KEEP=  0 

I F ( (KEY-LKEY) . LT . 1 )GO  TO  125 
LKEY=KEY 

WRI TE ( 6 , 1 ) I  TER , JMI N , TCOST , G 
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1  FORMAT (////IX, ' I TER= ' ,I4,5X,I4, ' COST  = ' ,F6.2,F12.8) 
WRITE (6, 4) 

4  FORMAT (//IX, ' U2= ' ) 

WRI TE (6,2) ( U2 ( I  ) ,1  =  1  ,50) 

2  FORMAT ( IX, 1 0F8 . 2  ) 

WRITE (6, 5) 

5  FORMAT (//IX,' GRAD=  '  ) 

WRITE (6, 3) (GRAD (I ) ,1=1 ,30) 

3  FORMAT ( IX, 1 0F8 . 2 ) 

WRITE (6, 6) 

6  FORMAT (//IX, ' TH= ' ) 

WRITE(6,3)  (TH(I  )  ,1  =  1  ,50) 

WRITE (6,8) 

WRITE (6, 3) (TL (I), 1=1, 50) 

8  FORMAT (//IX, ' TL= ' ) 

C . 

C  PART  THREE 

C  CONDITION  TO  STOP 

C . 

IF ( ABS (GM2-GMA2 ) /GM2 . LE . 0 . 1 )GO  TO  125 
STOP 


END 
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APPENDIX  III 


Program  to  simmulate  the  WATSUN  approach  and  the  suboptimal 

controller . 


C . 

C  THIS  PROGRAM  IS  USED  TO  SIMMULATE  THE 
C  SUBOPTIMAL  CONTROLLER  AND  THE  WATSUN  APPROACH 
C . 

c 

c 

c . 

C  TO  RUN  THE  PROGRAM: 

C  *USE  ADAPTIVE  CONTROLLER: 

C  CALL  READI N , WEATHR , AVERA , OPTSET , SYSANA 
C  *USE  WATSUN  APPROACH: 

C  CALL  READI N , WEATHR, SYSWAT 

C . 

COMMON  /INITL/THI , TLI ,U2I ,NDAYS 

COMMON  /I  FILE/I AMB , I  AVER, I  COST, I  DATA, IOPTI  , 

+  I  PARA, I  RAD 

COMMON  /WEATH/SLOPE , AZMTH , XLAT , RHO 

COMMON  /PARA/AB , AC , AAC , AH , ALPTAU , CC , CH , COPMAX , U 1  MAX , 
+  U2MAX , UL , TB , TMAX , ZETA , DT , FU3 
CALL  READI N 
CALL  WEATHR 
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CALL  AVERA 
10  CALL  OPTSET 
CALL  SYSANA 
C  CALL  SYSWAT 

STOP 
END 

C . 

BLOCK  DATA 

COMMON  /PARA/AB , AC , AAC , AH , ALPTAU , CC , CH , COPMAX , U 1 MAX , 
+  U2MAX , UL , TB , TMAX , ZETA , DT , FU3 
COMMON  /INITL/THI  , TLI  ,U2I  ,NDAYS 
COMMON  /I  FILE/IAMB, I  AVER, I  COST, I  DATA, IOPTI  , 

+  I  PARA, I  RAD 

COMMON  /WEATH/SLOPE , AZMTH , XLAT , RHO 
DATA  I AMB , I  AVER , I  COST ,  I  DATA , I OPTI  , I  PARA , 

+  I  RAD/2, 3, 7, 9, 4, 8, 1/ 

END 


C . 

C  THIS  SUBROUTINE  USED  TO  DETERMINE  THE  PURCHASED  ENERGY 
C  CONSUMED  BY  THE  SYSTEM,  USING  THE  SUBOPTIMAL  CONTROLLER. 
C . 

c 

SUBROUTINE  SYSANA 
DIMENSION  TH (200) , TL (200) 

REAL  MASS , MONEY, INSOL 
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COMMON  /PARA/AB , AC , AAC , AH , ALPTAU , CC , CH , COPMAX , 

+  U 1MAX , U2MAX , UL , TB , TMAX , ZETA , DT , FU3 
COMMON  /INITL/THI ,TLI ,U2I ,NDAYS 
COMMON  /I  FILE/IAMB , I  AVER , I  COST , I DATA , IOPTI  , 

+  I  PARA, I  RAD 

COMMON  /WEATH/SLOPE , AZMTH , XLAT , RHO 
C . 

10  FORMAT ( 1  OX , 2E 1 5 . 8 ) 

11  FORMAT ( 3  6X , E 1 5 . 8 ) 

C 1 2  0  FORMAT ( / ,  1 H  ,  ' D '  , I  4 ,  ’ COST=  ',E15.8,'  COP=',F6.2, 

C  +  /,’  REF.  COST= ' , E 1 5 . 8 , '  %SAVE= ' , F6 . 2 , / ) 

110  FORMAT( 1H  , ' S=’ ,E10.4, '  TA=’,F6.2,'  TH=’,F6.2, 

+  '  TL= ’  ,  F6 . 2 ,  ’  U2=  ’  , F5 . 2  ) 

Cl  60  FORMAT!/,  1H  ,  ’ ANNUAL  COST  =  ’  , E 1 5 . 8 ,  ’  ENERGY  COST= ’  , 
C  +  El  0.5,/,’  COL.  AREA= ’  ,F6 .  1  ,  ’  MASS  STORAGE^ ’  , E  1  0 . 5 ) 

C . 

REWIND  I  COST 
COLMAS  =  7  5 . 

I SAMP=  4 
5  CONTINUE 

C . 

C  SCALE  FACTOR  IS  1000 

C . 

AREA=AAC* 1 000 . 

MASS= ( CC* 1 000 . ) *2 . /4 . 18 
REWIND  I  DATA 


REWIND  IOPTI 
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C . 

C  INITIAL  CONDITONS 

C . 

I DAY= 1 
TH ( 1 ) =THI 
TL ( 1 ) =TLI 
U2=U2 I 
TOTAL=0 . 

RTOTAL=0 . 

SOLAR=  0 . 

I NSOL=  0 . 

C . 

C  SOLAR=SOLAR  COLLECTION 

C  I NSOL  =  I N SOL AT I  ON 

C  H POUT = OUTPUT  OF  HEAT-PUMP 

C  HPI N= INPUT  OF  HEAT-PUMP 

C  HPOUTD= OUTPUT  OF  HEAT-PUMP  IN  A  DAY 

C  HPI ND= INPUT  OF  HEAT-PUMP  IN  A  DAY 

C . 

HPOUTD=  0 . 

HPI ND=  0 . 

HPOUT=  0 . 

HPI N=  0 . 

I MARK=  0 
KOUNT= 1 
90  CONTINUE 
C  STOTAL=  0 . 
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C  DO  200  I TER= 1,10 

COST=  0 . 

RCOST=  0 . 

DO  150  1=2,24 
TH ( I )  =  0 . 

TL (I )  =  0 . 

150  CONTINUE 
THN  =  0  . 

C . 

READ ( I OPTI , 1 1 ) THN 
DO  100  1=1,24 
READ (I  DATA,  10)TA,S 

C . 

c 

C  I F ( I DAY . EQ . NDAYS ) GO  TO  29 

C  I F ( I DAY. NE. FOUNT) GO  TO  30 

C  29  .  WRITE (I  COST,  1 1 0 ) S , TA , TH ( I )  , TL ( I ) ,U2 
30  I F ( I MARK . NE . 1 ) GO  TO  40 

CALL  ONETAN ( S , TA , TH ( I )  , TL ( I )  , TH ( I  +  1  )  , 

+  TL ( I  +  1  )  , COST  1  , RCOST 1  , SOLAR  1  , I  SAMP ) 

U2  =  0  . 

HP=  0  . 

GO  TO  50 

C . 

I 

c 

40  CALL  TEMP ( S , TA , TH ( I ) , TL ( I ) , U2 , MARK , 

+  TH ( I  +  1 )  , TL ( I  +  1 )  , THN , COST  1  , HP , RCOST 1  , SOLAR  1  , I  SAMP ) 
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50  COST=COST+COST 1 

SOLAR=SOLAR  + SOLAR  1 
I NSOL= I NSOL+S 
RCOST=RCOST+RCOST 1 
HPOUTD=HPOUTD+HP 
HPI ND  =  HPI ND  +  U2  *U2MAX 

C  WRI TE ( I  COST ,  1  1 0 ) S , TA , TH ( I  )  , TL ( I )  , U2 

HPOUT=HPOUT+HP 
HPIN  =  HPI N  +  U2  *U2MAX 
S  =  0  . 

TA=  0  . 

I F ( TL ( I + 1 ) . NE . TH ( I + 1 ) ) GO  TO  179 
I F ( TL ( I  +  1 ) . GT . 2 4 .  ) GO  TO  100 
I MARK=  0 
GO  TO  100 

179  I F ( TL ( I  +  1 ) . LT .  ( TH ( I + 1 ) +  5 .  )  ) GO  TO  100 
AT= ( TL ( I + 1 ) +TH ( I + 1 ) )/2. 

TL (1+1 ) = AT 
TH (1+1 ) = AT 
I MARK= 1 
100  CONTINUE 

C  I F ( I DAY . EQ . NDAYS ) GO  TO  129 

C  I F ( I  DAY . NE . FOUNT ) GO  TO  130 

IF (HPIND.NE . 0 . ) COPD=HPOUTD/HPIND 
I F ( HPI ND . EQ . 0 . ) COPD=  2  0 . 

HPOUTD=  0 . 


HPI ND=  0 . 
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C  1  2  9 
1  1  9 
C 

C.  .  . 

130 

180 

C 

200 

c 

c 

400 

300 

c 

c 
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SAVE  = ( RCOST-COST ) *  1 00 . /RCOST 

WRI TE ( I  COST ,119)1  DAY , COST 

FORMAT ( / ,  1 H  ,  '  D  '  ,  I  4  ,  ' COST=  ’,E15.8,/) 

WRI TE ( I  COST ,12  0)1  DAY , COST , COPD , RCOST , SAVE 


TH ( 1 ) =TH ( 2  5 ) 

TL ( 1 ) =TL (25) 

I DAY=I DAY+ 1 

TOTAL=TOTAL+COST 

STOTAL  =  STOTAL  +  CO ST 

RTOTAL=RTOTAL+RCOST 

I F ( I DAY . GT . NDAYS ) GO  TO  300 

CONTINUE 

FOUNT =KOUNT+ 1 0 

WRI TE ( I  COST ,400  ) STOTAL 

FORMAT ( / , 1 H  , ' COST  OF  10  DAYS= ' , E 1 5 . 8 ) 

I F ( I  DAY . LE . NDAYS ) GO  TO  90 
CONTINUE 

WRITE( I  COST, 4  00  ) STOTAL 
WRITE (I  COST,  159) TOTAL 

FORMAT ( / ,  1 H  ,' SEASONAL  COST  ( KJ )  =  '  , E  1  5 . 8 ) 
COPHP=HPOUT/HPI N 

FORMAT (1H  ,’ TOTAL  ENERGY= ’  , E 1 0 . 5  ,  ’ COL .  AREA= '  , E 1 0 . 5 , 
+  ' MA  S  S  =  '  , E 1 0 . 5 , / ,  '  COP= '  ,  F4 . 2 ,  '  REF.  COST= '  , E 1 0 . 5 , 

+  ’  SYSCOP= '  , F6 . 4 ,  '  COLCOP= '  , F  6 . 2 , / ,  '  UL=',F6.2, 

+  ’  ALPTAU= ' , F4 . 2 , ’  U 1MAX= '  , F8 . 2 , '  AC=’,F4.1,/, 

+  ’  U2MAX= '  , F 6 . 2  ,  ' KW '  ,  '  ZETA= '  , 
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+  F4  .  1  ,  '  SLOPE= '  ,  F6 . 2 ,  '  DEG','  COLMAS=', 

+  F6 . 2 , / , '  NUMBER  OF  SAMPLES= ' ,13,'  PER  HOUR') 

S YSCOP= ( RTOTAL- TOTAL ) /RTOTAL 
COLCOP= SOLAR/INSOL 
ACCOEF=AC* 1000. 

UU 1 MAX=U 1 MAX* 1000. 

UU2MAX=U2MAX*  1  000  ./3600  . 

S S LOPE  =  SLOPE* 1 8  0 . / 3 .  1415927 

WRI TE ( I  COST ,16  1) TOTAL , AREA , MASS , COPHP , RTOTAL , SYSCOP , 
+  COLCOP , UL , ALPTAU , UU 1 MAX , ACCOEF , UU2MAX , ZETA , SSLOPE , 
+  COLMAS , I  SAMP 
RETURN 
END 


SUBROUTINE  TEMP ( S , TA , TH 1 , TL 1 , U2 , MARK , 

+  TH2 , TL2 , THN , COST  1  , HP, ROOST  1  , SOLAR  1  , I  SAMP) 

COMMON  /PARA/AB , AC , AAC , AH , ALPTAU , CC , CH , COPMAX , U 1 MAX , 
+  U2MAX , UL , TB , TMAX , ZETA , DT , FU3 
THF 1 =TH 1 
TLF 1 =TL 1 

DT= 1  . /FLOAT ( I  SAMP) 

TERM  1=0. 

COST  1=0. 

HP=  0  . 

RCOST 1=0. 


SOLAR  1 =0 . 


I TER= 1 


DELT=TB-TA 
20  TL 1 1 =TL 1 

100  I  F ( TH 1 . LE . TB ) U 1 =  0 . 

I F ( TH 1 . LE . TB ) GO  TO  111 
U 1 =AB*DELT/ ( U 1 MAX* ( TH 1 -TB ) *  ZETA ) 

1  1  1.  I  F  ( U  1  .GT.  1  .  )  U  1  =  1  . 

I F ( U 1 . LE . 0 . ) U 1  =  0 . 

QAUX=AB* (TB-TA)-UI *U1MAX* ( TH 1 -TB) *ZETA 
I F ( QAUX . LE . 0 . ) QAUX=  0 . 

UNUM=U1*U1MAX*ZETA* (TH1~TB)+AC* (TH1-TA) 

C . 

120  CONTINUE 

C . 

c . 

C  IF  TL1  . LT.  TA  USE  AMBIENT  AS 
C  SOURCE  FOR  HEAT  PUMP. 

C . 

I F ( TL 1 . LT . TA ) TL 1 1=TA 

UDEN=U2MAX* ( 1 . +2 . 5* ( 1 . - ( TH 1 -TL 1 1 ) /TMAX ) ) 
TERM2=AAC*FU3* ( S *ALPTAU~UL* ( TL 1 -TA ) ) 

I F ( TERM2 . LE . 0 .  ) TERM2  =  0 . 

I F ( TL 1 .GE.80. ) TERM2  =  0 . 

U2=0 . 5*U2+0 . 5*UNUM/UDEN* ( 1 . + ( THN-TH 1 ) ) 

101  CONTINUE 

COP= 1 .  +  ( COPMAX- 1 . ) * ( 1 . - ( TH 1 -TL 1 ) /4  5 . ) 

IF (COP.LT. ( 1 . +  AH/ ( ZETA*U 1 MAX ) ) )U2=0 . 


I F ( U2 . GT . 1 . ) U2= 1 . 

IF(U2.LE. 0. ) U2  =  0  . 

RH=UDEN*U2-UNUM 

I F ( TL 1  . LT . TA ) TERM  1 =AC* ( TL 1 -TA ) +TERM 1 
I F ( TL 1 . LT . TA ) GO  TO  107 

TERM 1 =  U2MAX*U2  *  2 . 5* ( 1 .-(TH1-TL1 1 ) /TMAX )  + 
+  AC* ( TL 1 -TA ) +TERM 1 

C . 

C  IF  TL 1  IS  LESS  THAN  TA ,  HEAT  PUMP  WON'T 
C  EXTRACT  HEAT  FROM  LOW-TEMPERATURE  TANK. 

C . 

107  RL=TERM2~TERM 1 

TH2=TH 1 +DT*RH/CH 

TL2=TL 1 +DT*RL/CC 

COST  1  =  ( QAUX+U2MAX*U2 ) *DT+COST 1 

HP  =  U2  *UDEN*DT+HP 

RCOST1=AB* (TB-TA) *DT+RCOST1 

SOLAR  1 =TERM2/AAC*DT  + SOLAR  1 

TL 1 =TL2 

TH 1 =TH2 

I TER= I TER+ 1 

I F ( I  TER . LE . I  SAMP ) GO  TO  2  0 
TH 1 =THF 1 
TL 1 =TLF 1 
RETURN 


END 
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C  THIS  IS  THE  SUBROUTINE  TO  READ 
C  PARAMETER  DATA  OF  THE  SYSTEM. 

SUBROUTINE  READIN 

COMMON  /PARA/AB , AC , AAC , AH , ALPTAU , CC , 

+  CH , COPMAX , U 1 MAX , U2MAX , UL , TB , TMAX , ZETA , DT , FU3 
COMMON  /INITL/THI ,TLI ,U2I ,NDAYS 
COMMON  /I F I LE/I AMB , I  AVER , I  COST , I  DATA , I OPTI  , 

+  I  PARA, I  RAD 

COMMON  /WEATH/SLOPE , A ZMTH , XLAT , RHO 
10  FORMAT(//, 3E1 0.3) 

20  FORMAT (/,8E1 0.3) 

30  FORMAT (/,4E 10. 3) 

40  FORMAT!/, 3E1 0 . 3 ) 

50  FORMAT!/, El  0.3) 

60  FORMAT !/,3E10.3,I5) 

C 

c . 

REWIND  I  PARA 

C . 

C  STORAGE  DATA 

C . 

READ ( I  PARA ,  1 0 ) AC , AH , U 1 MAX 
C 

c . 

C  COLLECTOR  DATA 

C . 


READ ( I  PARA, 20) AAC , ALPTAU , FU3 , UL , SLOPE, XLAT , 
+  AZMTH  , RHO 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BUILDING  DATA 


READ ( I  PARA , 3  0 ) AB , CB , TB , ZETA 


HEAT-PUMP  DATA 


READ ( I  PARA ,40) COPMAX , TMAX , U2MAX 


SAMPLING  TIME  (STEP  TIME) 
READ ( I  PARA , 50 ) DT 


INITIAL  CONDITIONS 


READ( I  PARA, 60 )THI ,TLI  ,U2I  , NDAYS 


SCALED  DATA (TO  MINIMIZE  ROUND-OFF  ERROR) 


C  IN  THE  FOLLOWING  ,  THE  STORAGE  IS  ASSUMED  TO  BE 
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C  7 5KG  PER  COLLECTOR  SQUARE  METER 

C . 

CC  =  4 . 1 8  *AAC*  7 5  . 

CH  =  CC 

C . 

AC=AC/1 000 . 

AH=AH/1 000 . 

U 1 MAX=U 1 MAX/ 1000. 

AAC=AAC/1 000 . 

AB=AB/1 000 . 

U2MAX=U2MAX*3600 . /I  00  0 . 
CC=CC/1 000 . 

CH=CH/1 000 . 

RETURN 

END 


C . 

C  SUBROUTINE  OPTSET  USED  TO  DETERMINE 
C  ALL  OPTIMUM  SETS  OF  TH , TL  BASED  UPON 
C  THE  AVERAGE  VALUE  OF  AMBIENT  TEMPERATURE 
C  AND  RADIATION. 

C . 

c 


SUBROUTINE  OPTSET 
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COMMON  /PARA/AB , AC , AAC , AH , ALPTAU , CC , CH , COPMAX , U 1 MAX , 
+  U2MAX , UL , TB , TMAX , ZETA , DT , FU3 
COMMON  /I  FILE/IAMB , I  AVER , I  COST , I  DATA , IOPTI  , 

+  I  PARA, I  RAD 

COMMON  /INITL/THI , TLI ,U2I ,NDAYS 
C 

10  FORMAT (7X, 2E15.8) 

15  FORMAT (1H  , ' S=  ' ,E15.8,?  TA=  ' ,F6.2,’  TH0  =  ’ ,E15.8) 

C 

c . 

REWIND  I  AVER 
REWIND  IOPTI 

C . 

DO  20  1=1 , NDAYS 

READ (I  AVER,  1 0 ) TA , S 

TH0  =  TB  +  AB* ( TB-TA ) / ( U 1 MAX*  ZETA ) 

WRITE (IOPTI , 1 5 ) S , TA , THO 
S  =  0  . 

TA=  0  . 

20  CONTINUE 
RETURN 
END 

C************************************************ 
SUBROUTINE  WEATHR 

C  . TO  PROCESS  METEOROLOGICAL  FILES 


U*  ^  'l'  X-  >X*  •l'  ^  *X-  ^  *1*  «I>  •X*  «X-  X#  X>  X«  X*  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X 

T*  *T*  *T-  'T  *T'  T'  V  *T'  T>  «T*  ^  T*  «T*  *T'  T  'T'  *T*  T  *TS  *T*  V  *T*  V  T'  ^  ^  4*  T*  T»  T*  ^  ^  ^  ^  ^  ^  ^  ^ 

INTEGER  RADTN (24) ,  I T  (  2  4  ) 

COMMON  /I  FI LE/I AMB , I  AVER , I  COST , I  DATA , I  OPT I , 
+  I  PARA, I  RAD 

COMMON  /WEATH/SLOPE , AZMTH , XLAT , RHO 
COMMON  /INITL/THI ,TLI ,U2I ,NDAYS 
FORMAT ( 17X, 24(16, IX)) 

FORMAT ( 17X, 24(16, IX)) 


22 
25 
2  1 
3  1 


FORMAT ( 1 5X, I  2 ) 

FORMAT ( 3X , ’ D ’  ,  I  3  ,  ’  H  ’  , I  2 , 2E 1 5 . 8 ) 


100  FORMAT ( 30X , 4  OX ) 

SINSLP=SIN( SLOPE ) 

COSSLP=COS( SLOPE ) 

S I NAZM=  S I N ( AZMTH ) 

COSAZM=COS (AZMTH) 

SINLAT=SIN (XLAT ) 

COSLAT=COS (XLAT) 

TERM  1  =  SI NLAT* COS SLP- COSLAT *SI NSLP*COSAZM 
TERM2 = COSLAT* COS SLP+ SI NLAT*SI NSLP*COSAZM 
TERM3=SI NSLP*SI NAZM 
RD1=( 1 . +COSSLP ) /4 . 

RD2= ( 1 . -COSSLP ) / 2 . 


REWIND  IAMB 


REWIND  I  RAD 


REWIND  I DATA 


C  TO  FORWARD  ONE  RECORD , COMMENT ,  OF  I  RAD  AND 


C  IAMB  DATA  FILES 


READ (IAMB, 100) 

READ ( I  RAD,  100) 

C  . 

DO  11  I DAY= 1 , NDAYS 

C .  30  SHOULD  BE  CHANGED  TO  365  OR  366  (DAYS) 

DAY=I DAY  +  27  3 

C .  OCTOBER  1ST  IS  THE  274TH  DAY  OF  THE  YEAR 

I F ( DAY . GT . 36  5  . ) DAY  =  DAY~  3  6  5 . 

41  READ ( I AMB ,21) KCRI T 

IF (KCRIT.EQ. 78 )GO  TO  42 
KCRI T=  0 
GO  TO  41 

42  BACKSPACE  IAMB 
READ ( IAMB , 22 ) IT 

4  3  READ ( I  RAD ,21) KCRI  T 

IF (KCRIT.EQ. 6 1 ) GO  TO  44 
GO  TO  43 

44  BACKSPACE  I  RAD 

READ ( I  RAD ,22) RADTN 
RHO=RHO*RD2 

SC  =  4 8 7  1  .  0*(  1 .+0.3 3 *COS ( 1  .72  1  42E-02*  DAY) ) 
DECL=0 . 40928*SIN ( 1 . 72 1 42E-02* ( 284 . +DAY) ) 

S I NDEC=SI N ( DECL ) 

COSDEC=COS (DECL) 

TERM4=S I NDEC*TERM 1 


TERM5=COSDEC*TERM2 


TERM6=COSDEC*TERM3 


TERM7  =  SI NDEC*S  I NLAT 
TERM 8= COS DEC* CO SLAT 

C . 

DO  1 0  I HR= 1 ,24 
TA= I T ( I  HR ) *  0 .  1 
H  =  RADTN ( I  HR ) 

C . 

I F ( H . LE . 0 . ) GO  TO  8 
HRANG=0 . 26 1 8* ( 1 1 . 5- 1  HR) 

S I NHR  =  S I N ( HRANG ) 

COSHR=COS ( HRANG ) 

C . 

C0STT  =  TERM4  +  C0SHR*TERM5  +  S I NHR*TERM6 
I F ( COSTT . LE . 0. ) GO  TO  8 

C . 

COSTZ=TERM7+COSHR*TERM8 
IF (COSTZ .LE. 0 . ) GO  TO  8 

C . 

RB= COSTT/ COSTZ 

IF( (RB.LE.O. ) .OR. (RB.GT.5. ) )RB=0. 
HEX= SC* COSTZ 

I F ( HEX . LT . H ) HEX=H 
XKT=H/HEX 

c . 

IF (XKT.LT. 0 . 35 )GO  TO  5 
IF (XKT.LT. 0 . 75) GO  TO  6 
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HD=H*  0 . 1769 
GO  TO  7 

5  HD  =  H* (  1 .0-0.248857  *XKT ) 

GO  TO  7 

6  HD=H* ( 1 .55699-1 .8401 3*XKT ) 

7  CONTINUE 


HB=H-HD 

C  HBT=(0.5*HD+HB)*RB 

C  HDT=HD*RD 1 +H*RHO 

C  HT=HBT+HDT 

HT  =  RB*HB  +  2 . *RD 1 *HD  +  RHO*RD2  *H 
GO  TO  9 

8  HT=  0 . 

9  WRITE  ( I  DATA  ,31)1  DAY  ,  I  HR  ,  TA  ,  HT 

10  CONTINUE 

11  CONTINUE 
END  FILE  I  DATA 

RETURN 

END 


C . 

C  THIS  PROGRAM  USED  TO  DETERMINE  THE 
C  AVERAGE  VALUES  OF  AMBIENT  TEMPERATURE 


C  AND  RADIATION. 


1  37 


C . 

SUBROUTINE  AVERA 

DIMENSION  T( 24 ) , S ( 24 ) ,TAV( 200 ) , SAV( 200 ) 
COMMON  /I  FILE/IAMB, I  AVER, I  COST, I  DATA, I OPTI , 
+  I  PARA, I  RAD 

COMMON  /INITL/THI , TLI ,U2I ,NDAYS 

c . 

REWIND  I  DATA 
REWIND  I  AVER 

C . 

DO  10  1  =  1, NDAY  S 
SAV ( I ) =0 . 

TAV ( I ) =0 . 

DO  1 1  J= 1 , 24 

READ ( I  DATA ,20)T(J),S(J) 

20  FORMAT ( 1  OX , 2E 1 5 . 8 ) 

TAV (I ) =  TAV ( I )+T(j) 

SAV ( I )=SAV( I )  +S ( J ) 

T ( J ) =0 . 

S  ( J  )  =  0  . 

11  CONTINUE 

TAV ( I ) =TAV ( I )/24 . 

SAV (I ) =  SAV ( I )/24. 

WRI TE ( I  AVER ,30)1, TAV ( I )  , SAV ( I ) 

3  0  FORMAT ( 3X , ' D '  , I  3 , 2E  1  5 . 8 ) 

10  CONTINUE 


RETURN 


138 


END 


C  THIS  PROGRAM  USED  TO  SIMMULATE  THE 
C  APPROACH  ADOPTED  BY  WATSUN  PROGRAM. 

C  HERE  HEAT  PUMP  IS  USED  IN  ON-OFF 
C  MODE  ONLY. 

SUBROUTINE  SYSWAT 
DIMENSION  TH ( 200 ) , TL ( 200 ) 

COMMON  /PARA/AB , AC , AAC , AH , ALPTAU , CC , CH , COPMAX , U 1 MAX , 
+  U2MAX , UL , TB , TMAX , ZETA , DT , FU3 
COMMON  /INITL/THI ,TLI ,U2I ,NDAYS 
COMMON  /I  FILE/IAMB , I  AVER, I  COST, I  DATA , IOPTI  , 

+  I  PARA , I  RAD 

WRITE ( I  COST,  12  1) 


REWIND  I  DATA 

C . 

C . 

C  INITIAL  CONDITIONS 

C . 

I SAMP=4 
I DAY= 1 
TH ( 1 ) =THI 
TL ( 1 ) =TLI 


HPOUT=0 . 
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HPI N=  0 . 

HPI ND=0 . 

HPOUTD=  0 . 

TOTAL=0 . 

KOUNT= 1 
90  CONTINUE 
STOTAL=  0 . 

DO  300  I TER= 1,10 
COST=  0 . 

DO  150  I T=  2 , 2 5 
TH ( I T ) =  0 . 

150  TL ( I T ) =  0 . 

DO  100  1=1,24 
READ (I  DATA,  1 0 ) TA , S 
10  FORMAT ( 1  OX , 2E 1 5 . 8 ) 

C  IF ( I  DAY .EQ.NDAYS )GO  TO  29 

C  IF (IDAY.NE. FOUNT ) GO  TO  30 

C2  9  WRITE (I  COST,  1 1 0 ) S , TA , TH ( I )  , TL ( I 
30  CALL  TEMPI ( S , TA , TH ( I ) , TL ( I ) , 

+  TH ( I  +  1 )  , TL ( I  +  1 )  , THN , COST  1  , HP , U2 
S  =  0  . 

TA=  0  . 

COST=COST+COST 1 
HPOUT=HPOUT+HP*CC/ 1000. 
HPIN=HPIN+U2*U2MAX*CC/ 1000. 
HPIND=HPIND+U2*U2MAX*CC/1 000 . 
HPOUTD=HPOUTD+HP*CC/1 000 . 


,U2 

I  SAMP) 


100  CONTINUE 


C . 

121  FORMAT ( // ,  1 H  ,  1  OX WATSUN  APPROACH ',/ ) 
C  I F ( I  DAY . EQ . NDAYS ) GO  TO  129 

C  I F ( I DAY . NE . FOUNT ) GO  TO  130 

IF (HPIND.NE. 0 . ) COPD=HPOUTD/HPI ND 
IF (HPIND.EQ. 0 . ) COPD=  2  0 . 

HPOUTD=  0 . 

HPI ND=  0 . 

C 1 2  9  WRI TE ( I  COST ,119)1  DAY , COST 
119  FORMAT ( / ,  1 H  , ’ D ’  , I  4 ,  ' COST=  ’,E15.8) 

C 1 2  9  WRITE (I  COST,  1 2 0 ) I  DAY , COST , COPD 

C . 

130  TH ( 1 ) =TH (25) 

TL( 1 ) =TL (25) 

I DAY= I DAY  + 1 
TOTAL = TOTAL + COST 
S TOTAL = S TOTAL + CO ST 
IF ( I  DAY .GT. NDAYS )GO  TO  3 1 0 
300  CONTINUE 
C  KOUNT=KOUNT+ 1 0 

WRITE ( I  COST ,400  ) STOTAL 
400  FORMAT ( / ,  1 H  , ’ COST  OF  10  DAYS= ’  , E 1 5 . 8 ) 
IF (IDAY.LE. NDAYS )GO  TO  90 
310  COP =H POUT /HP I N 

WRITE( I  COST, 40  0  ) STOTAL 
C  WRITE ( ICOST, 20 1 ) TOTAL 


201  FORMAT ( / , 1 H  , ’ SEASONAL  COST  ( K J ) = ' , E 1 5 . 8 ) 

WRITE ( I  COST, 200  ) TOTAL, COP 

200  FORMAT ( / , 1 H  ,  '  TOTAL  COST=  ',E15.8,'  COP=’,F6.2/) 

120  FORMAT ( / , 1 H  , '  D  ’ ,13, IX, ’ COST=  \E15.8,’  COP=’,F6.2/) 
110  FORMAT ( 1 H  ,  ’ S= ’  ,E  1  0 . 4  ,  ’  TA=’,F6.2,’  TH=',F6.2, 

+  ’  TL= ’ , F6 . 2 , ’  U2= ’ , F6 . 2 ) 

RETURN 

END 

SUBROUTINE  TEMP  1 ( S , TA , TH 1  , TL 1  , 

+  TH2 , TL2 , THN , COST  1 , HP , U2 , I  SAMP ) 

COMMON  /PARA/AB , AC , AAC , AH , ALPTAU , CC , CH , COPMAX , U 1 MAX , 
+  U2MAX , UL , TB , TMAX , ZETA , DT , FU3 
THF 1 =TH 1 
TLF 1 =TL 1 

DT= 1  . /FLOAT ( I  SAMP) 

COST  1=0. 

HP=  0  . 

TERM  1  =  0. 

DELT=TB-TA 
I TER= 1 

20  TL 1 1 =TL 1 

IF (TL1 .LT. 20 . ) GO  TO  100 

C . 

C  IF  LOW-TEMPERATURE  TANK  COULD  PROVIDE 
C  HEAT  TO  BUILDING,  IT'S  WELCOME. 
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C . 

U1P=AB*DELT/(U1MAX* (TL1-TB) ) 

I F ( U 1 P . GT . 1 . ) U 1 P= 1 . 

I F ( U 1 P . LT . 0 . )U1P=0. 

QAUX=AB* (TB-TA)~U1P*U1MAX* (TL1-TB) 

I F ( QAUX . LE . 0 . ) QAUX=  0 . 

IF (QAUX.NE. 0 .  )G0  TO  115 

C . 

C  IN  CASE  LOW-TEMPERATURE  TANK  CANNOT  PROVIDE 
C  ENOUGH  HEAT,  HIGH-TEMPERATURE  TANK  WOULD 
C  COME  IN. 

C . 

UNUM=AC* ( TH 1 -TA ) 

TERM  1 =U 1 P*U 1 MAX* (TL1-TB) 

GO  TO  120 

115  U1=QAUX/(U1MAX* (TH1-TB) ) 

I F ( U 1 .GE. 1 . ) U 1 = 1 . 

IF (U 1 .LT. 0 . )U 1 =0 . 

QAUX = QAUX -U 1  *U  1 MAX* (TH1-TB) 

I F ( QAUX . LE . 0 . ) QAUX= 0 . 

UNUM=U 1 *U 1MAX* ( TH 1 -TB ) +AC* ( TH 1 -TA ) 

TERM  1 =U 1 P*U 1 MAX* (TL1-TB) 

GO  TO  120 

100  U 1 =  AB  *DELT/ ( U 1 MAX* ( TH 1 -TB )  ) 

IF ( U 1 .GT. 1 . ) U 1 = 1 . 

IF (U1 . LT . 0 . ) U 1 =0 . 

QAUX=AB*(TB-TA)-U1*U1MAX*(TH1-TB) 


I F ( QAUX . LE . 0 . ) QAUX= 0 . 

UNUM=U 1 *U 1 MAX* ( TH 1 -TB ) +AC* ( TH 1 -TA ) 

120  CONTINUE 

C . 

C  IF  TL 1  .LT.  TA  USE  AMBIENT  AS 
C  SOURCE  FOR  HEAT  PUMP. 

C . 

I F ( TL 1 . LT . TA ) TL 1 1=TA 

UDEN  =  U2MAX* (  1  . +2 . 5* (  1 . - ( TH 1 -TL 1  1 )/TMAX) ) 
TERM2 = AAC*FU3 * ( S*ALPTAU~UL* ( TL 1 -TA ) ) 

I F ( TERM2 . LE . 0 . )TERM2=0. 

C . 

C  SOLAR  ENERGY  DUMPED  IF  TL1>20. 

C . 

IF (TL1 .GT. 2  0 . ) TERM2  =  0 . 

C . 

C  OTHERWISE,  TURN  ON  HEAT-PUMP, 

C  ELECTRICAL  INPUT  TO  THE  HEAT  PUMP  IS 
C  APPROXIMATED  BY  THE  FOLLOWING  EXPRESSION 

C . 

U2= ( 2 . + ( TL 1 -TH 1 +44 . )/25. )*3.6/U2MAX 
IF (U2 .GT. 1 . ) U2= 1 . 

IF (U2 .LT. 0 . )U2=0 . 

RH=UDEN*U2-UNUM 

I F ( TL 1 . LT . TA ) TERM 1  =AC* ( TL 1 -TA ) +TERM 1 
I F ( TL 1 . LT . TA ) GO  TO  107 

TERM  1 =U2MAX*U2* 2 . 5* ( 1 .- (TH1-TL1 1 ) /TMAX )  + 


144 


+  AC* ( TL 1 -TA ) +  TERM 1 

C . 

C  IF  TL 1  IS  LESS  THAN  TA ,  HEAT  PUMP  WON'T 
C  EXTRACT  HEAT  FROM  LOW-TEMPERATURE  TANK. 

C . 

107  RL=TERM2-TERM1 

TH2=TH 1 +DT*RH/CH 

TL2=TL 1 +DT*RL/CC 

COST  1  =  ( QAUX+U2MAX*U2 ) *DT  +  COST 1 

HP  =  U2  *UDEN*DT+HP 

TL 1 =TL2 

TH 1 =TH2 

I TER= I TER+ 1 

IF ( ITER. LE. I  SAMP) GO  TO  2  0 

TH 1 =THF 1 

TL 1 =TLF 1 

RETURN 

END 


SUBROUTINE  ONETAN ( S , TA , TH 1  , TL 1  , TH2 , TL2 , COST  1 , RCOST 1  , 
+  SOLAR  1  , I  SAMP) 

C . 

C  THIS  SUBROUTINE  USED  IN  THE  CASE  WHERE  TWO  TANKS  ARE 
C  TO  BECOME  ONE-TANK  SYSTEM, I. E.  S  HIGH  AND  TA  LOW 
C . 


145 


COMMON  /PARA/AB , AC , AAC , AH , ALPTAU , CC , CH , COPMAX , U 1 MAX , 
+  U2MAX , UL , TB , TMAX , ZETA , DT , FU3 
THF=TH 1 

DT= 1  ./FLOAT (I  SAMP) 

I TER= 1 
SOLAR  1=0. 

RCOST 1=0. 

COST  1=0. 

DELT=TB-TA 
TT 1 =TH 1 
10  CONTINUE 

U 1 =AB*DELT/ ( U 1 MAX* (TT1-TB) ) 

I F ( U 1 .GT. 1 . ) U 1 = 1 . 

I F ( U 1 . LT . 0 . ) U 1 =0 . 

QAUX=AB*  (TB-TA)  -UUU1MAX*  (TT1-TB) 

IF (QAUX.LE. 0 .  ) QAUX=  0 . 

UNUM=U 1 *U 1 MAX* ( TT 1 -TB ) +  AH* ( TT 1 -TA ) + 

+  AC* ( TT 1 -TA ) 

TERM=AAC*FU3* ( S*ALPTAU~UL* (TT1 -TA) ) 

I F ( TERM . LE . 0 . ) TERM= 0 . 

IF (TT1  .GE. 90 . ) TERM=  0 . 

RTT=TERM-UNUM 

TT2 =TT 1 +DT*RTT/ (CC+CH) 

COST  1 =QAUX*DT  +  COST 1 

RCOST 1=AB* (TB-TA) *DT+RCOST1 

SOLAR  1 = TERM/A AC  *DT  + SOLAR  1 


I TER= I TER+ 1 
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TT 1 =TT2 

IF ( ITER. LE. I  SAMP) GO  TO  10 

TH 1 =THF 

TL 1 =THF 

TH2=TT2 

TL2=TT2 

RETURN 

END 


appendix  IV 


Properties  of  the  system  eigenvalues 


The  dynamic  equations  for  the  solar-assisted  heat  pump  system  are 


(equations  (2.3-la)  and  (2.3-lb): 


T  -T 

C  T  =  -£u  (T  -T  )  -  a,  (T  -T  )  +  u_[l  +  (COP  -1)  (1  -  - - -)  ] 

hh  l  hb  hha  2  max  T 


max 


T  -T 

C  T  =  -u  r(COP  -1) (1  -  - - -)]  -  a  (T  -T  )  + 

c  c  2  max  T  c  c  a 

max 


AcF(u3)[S  -  UL(Tc  -  V1 


or : 


(COP  -1)  u  (COP  -1) T 

_  r  _  max  ,  ^  ,  2  max  c  , 

ChTh  -  [-5ul  -  \  -  U2  T - ]Th  + - 5 -  + 


max 


max 


5u  T  +  a  T  +  u  COP 
lb  ha  2  max 


(IV. 1) 


C  T 
c  c 


u 

2 


(COP  -1) 
max 


T 

max 


Th  + 


[-u 


(COP  -1) 
max 

T 


max 


a 

c 


A  F (u  )U  ] T 
c  3  L  c 


-u0 (COP  -1)  +  a  T  +  A  F (u  ) (S  +  U  T  ) 


2 


max 


c  a 


L  ay 


Let  X 


~ 

—  f 

X1 

T,  ' 

!  l 

h 

T 

i  2 

c 

L  - 

—  — 

(IV. 2) 
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"?V  ah  '  u2 


4  - 


u0 (COP  -1) 
2  max 

T  C 
max  c 


(COP  -1) 
max 


max 


u_(COP  -1) 
2 _ max 

T  CL 
max  h 


U2(C°Pmax  -  a  -  A  F(u  )U 

- - -  c  c  3  L 

max 


(uT  +  a  T  +  uoC0P 
Ip  ha  2  max 


r  = 


-u_ (COP  -1)  +  a  T  +  A  F (u  ) (S  +  U  T  ) 
2 _ max  c  a _ c  3  La 

C 


Equations  (IV. 1)  and  (IV. 2)  become: 


x  =  £  x.  _L 


(IV. 3) 


In  the  following,  we  will  prove  that  the  eigenvalues  e.g.  (IV. 3)  are 
always  nonpositive. 


Proof :  According  to  the  Gerschgorin  theorem  [l],if  we  designate 

A  A„  as  the  eigenvalues  of  the  system  described  by  (IV. 3),  we  obtain: 

_L ,  Z 


(COP  -1) 

r  >-  max 

(-5ul  '  ah  -  U2 - T - 1 


max 


u  (COP  -1) 
2 _ max 

T  CL 
max  h 


(IV. 4) 


149 


a  nd , 


X2 


'U2<'C0Pmax  1)  -  a  -  A  F(u  )U 

- ^ -  c  c  3  L 

max 


u0  (COP  -1) 
2  max 

T  C 
max  c 


c  (COP  -1) 

u.  -  a.  -  u.  max 
I  n  2- 


Let  - 


T 
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,-u  (COP  -1)  ,  \ TT 

2  max  -a  -  A  F(u„)UT 
-  c  c  3  L 


f2  ^  A2  - 


T 


max 


If  f,  >  0  since  COP  -  1  >  0  and  un  >  0,  (IV. 4)  becomes 
1  -  max  2 


- 


a,  -  u  (COP  -1) 
h  2  max 


b  - 


max 


.  u  (COP  -1) 
<  __2 _ max 

T  C 
max  h 


or  X  <  -Cui  +  ah) 

Ch 

Hence,  is  nonpositive. 

If  f  <  0,  (IV. 6)  becomes: 


(IV. 5) 


(IV. 6) 


(IV. 7) 


(IV. 8) 
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-Cu,  - 


a  -  u  (COP  -1) 
h  2  max 


fi  -  h  - 


max 


<  0 


or 


+  a,  +  u  (COP  -1) 
1  h  2  max 


h  i 


max 


Therefore,  A^  is  always  nonpositive. 

Similarly,  we  can  also  prove  that 


(a  +  A  F (u  )U  ) 
^2  -  ~  c  c  3  L 


(IV. 9) 


(IV. 10) 


or 


+  A  F (u  )U 
c  3  L 


C 

c 


+  u  (COP  -1) 
2 _ max 

T 

max 


(IV. 11) 


If  and  a^  are  different  from  zero,  from  (IV. 8),  (IV. 9),  (IV. 10), 
(IV. 11),  we  find  that  the  eigenvalues  of  the  system  are  always  negative 
no  matter  what  the  value  of  u^,  u?,  F(u  )  are. 

Let  u^  =  F(u^)  -  1,  matrix  <£  reduces  to: 


^  =  A  +  B  u2 

where  A  and  are  coefficient  matrices  of  system  (3.1-10).  Hence,  with 

nonzero  values  of  a,  and  a  ,  the  eigenvalues  of  the  system  described  by 

h  c 

x  =  <£  _x  =  (A  +  J3  u^)  _x 

are  always  negative  regardless  of  the  values  of  u?. 
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APPENDIX  V 


Deadbeat  Control  of  Bilinear  Plant 

In  the  following,  the  controllability  and  optionality  problems 
of  a  general  discrete-time  bilinear  system  are  investigated.  Although, 
the  results  look  promising  from  the  viewpoint  of  the  suboptimality,  it 
is  imperative  that  more  work  needs  be  done  since  there  are  still  unsolved 
problems  with  this  approach  (for  example,  there  are  problems  of  imple¬ 
menting  the  controller  as  well  as  of  setting  up  the  state-transition 
matrix,  etc.) 

For  a  continuous  bilinear  plant,  its  behaviours  can  also  be  examined 
in  terms  of  this  approach,  since,  as  we  will  show,  the  discrete-time 
equivalence  of  a  continuous  bilinear  plant  is  also  discrete-time  bilinear. 
Assuming  that,  the  general  continuous  bilinear  system  equations  can  be 
written  as  follows: 


X(t) 


N  N 

u .  v . 

1  1 

(H  +  £  u.E  +  £  v.  E  ) x ( t )  +Du+Dw+D 

—  .  _  1— U  .  .  1.  —V  .  —  — 'V-  -^W—  — o 

1=1  1  1=1  1 
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—  -  — v  —  -v?  —  — o 


(V.l) 
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l 
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N 
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l 


u  E  + 
l  — u . 
l 


L 

i=l 


v.  E 
i  ~~v . 
l 
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H,  E  ,  E  ,  D  ,  D  ,  D  =  constant  matrices 
—  —  i  .  —v .  —v  — w  — o 
1  l 


tj,  v,  w  =  control  vectors 

The  homogeneous  solutions  to  the  above  differential  equation  is 


F(t-t  ) 

25^(0  =  e  x(t:  ) 


where,  by  definitions,  the  exponential  matrix  is  equal  to 


F(t-t  )  0/  .2 

e  °  =  I  +  F (t-t  )  +  F2(t  to)  + 

=  —  o  — 


2  ! 


or  by  linear  approximation,  we  have 


F(t-t  ) 

e  -  I  +  F (t-t  ) 

—  o 


Hence , 


jc  (t)  =  [  I  +  F ( t- 1  )  ]  x  (t  ) 

— h  —  o  —  o 


(V .  2) 


The  particular  solution  to  (V.l)  is  obtained  by  using  the  familiar 
method  , i . G . variation  of  parameters.  Assuming  that,  the  solution  is 
in  the  form 


F(t-t  ) 

x  (t)  =  e  z(t) 

— p  - 


(V .  3) 


Substituting  (V.3)  into  (V.l),  we  obtain 


F(t-t  )  F(t-t  ) 

Fe  z(t)+e  z(t)=Fx(t)+D  v+D  w+D 

—  —  —  - ^p  —v  —  w  —  — o 


F(t-t  )  F(t-t  )  F(t-t  ) 

Fe  °  z(t)  +  e  z(t)  =  Fe  °  z(t)  +D  v+D  w+D 

—  — v  '  —  —  —  —v  —  — w  —  — o 


or 
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This  implies 


F(t-t  ) 

e  °  z(t)  =  D  v+D  w  +  D 

—V  —  “W  —  — o 


so 

z(t) 


-F(t-t  ) 

e  °  (D  v  +  D  w  +  D) 

—v  —  — w  —  — o 


z(t) 


t  to) 

=  /  e  (Dv+Dw+D)d£ 

—v  —  — w  —  — o 

t 

o 


(V . 4a) 


Suppose,  vector 


(D  v+D  w+D)  does  not  depend  on 
—v  —  w  —  — o 


(V.4a)  can  be 


written 


t  —(+0^ 

z(t)  -f  e  d£  (D  v  +  D  w  +  D  ) 

—  ~^v  —  — w  —  — o 

t 

o 


(V . 4b) 


Thus,  from  (V.4b)  and  (V.3) 


x  (t) 
“ P 


F(t-to) 
e  z  (t) 


t  -FU-t  ) 
/  e 
t 


dr(D  v  +  D  w  +  D) 
^  —iT  —  — w  —  — o 


(V .  5 ) 


The  complete  solution  to  (V.l)  is 


x(f)  =  xh(t)  +  Xp(t) 


Since,  we  wish  to  use  this  solution  over  one  sample  period,  let 

t  =  nT  +  T  and  t  =  nT,  we  have: 

o 
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x(nT  +  T)  =  x,  (nT  +  T)  +  x  (nT  +  T) 
—  — h 


nT+T  F (nT+T-£) 

=  (I  +  FT)  x(nT)  +  /  e  d? (°v  X  + 

nT 


In  equation  (V.6),  we  assume  that  (using  zero-order  hold) 

u(0  =  u_(nT)  for  nT  <  £  <  nT  +  T 

v(c)  =  v(nT)  nT  <  L,  <  nT  +  T 

w(^)  =  w(nT)  nT  <  £  <  nT  +  T 


nT+T  F(nT+T-C) 


Let  F  =  /  e 

1  nT 


dC  with  shortened  notation,  equat 


becomes  x.(n+l)  =  (I_  +  _FT)  }£(n)  +  Z.n  (D  v(n)  +  D  w(n)  +  D 

1  V  w 


N 


N 


u 


v 


or,  x (n+1)  =  (I  +  HT  +  E  u. (n)  E  T  +  E  v.(n)  E  T) 

~  —  —  .  1  — U.  .  -i  i  “V. 

i=l  i  i=l  i 


F.  D  v (n)  +  F,  D  w(n)  +  F,  D 

—1  —v  —  —1  — w  —  —1  — o 


If  we  let 


A  =  I  +  HT 


B  =  E  T 
— u  .  — u  . 

l  l 


B  =  E  T 
— - v .  —v . 

l  l 


C  =  F  D 
— v  — 1  ~^v 


C  =  F  D 

-~w  —1  ~“W 


C  =  F  D 
— o  —1  — o 


D  w  +  D  ) 

— W  —  — o 

(V.6) 


ion  (V.6) 

)  (V .  7) 
o 

x(n)  + 

(V .  8) 

(V .  9) 
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and  change  n  by  k,  we  obtain 


N  N 

u  v 

x(k+l)  =  (A  +  Z  u. (k)  B  +  Z  v. (k)  B  )  x(k)  + 
—  —  x  ~u.  .  ,  1  —v.  — 


C  w(k)  +  C  v(k)  +  C  (V.10) 

— W  —  — V  —  — o 


The  above  is  a  general  discrete- time  belinear  system  equation.  Thus, 
it  is  possible  to  approximate  to  the  first  order  a  continuous  bilinear 
plant  by  an  equivalent  discrete-time  bilinear  one. 

The  general  discrete- time  plant  equation  is  rewritten  as  follows 

[1]: 


N  N 

u  v 

x(k+l)  =  [A  +  Z  u.(k)  B  +  Z  v . (k) 
—  —  .  ,  l  — u .  .  ,  i 


B  ]  x (k)  +  C  v (k)  + 


C 


w(k)  +  C 
—  —  o 


(V.ll) 


where , 

x(k)  =  state  vector  at  time  (k) 
x(k+l)  =  state  vector  at  time  (k+1) 

(k)  =  multiplicative  control  input  (i)  at  time  (k) 

a  - 1,  ....  y 

v  (k)  =  multiplicative  and  additive  control  input  (i)  at  time  (k) 
i 

(i  =  1,  Nv) 

w(k)  =  additive  control  input  vector  at  time  (k) 

C  =  vector  of  constant  inputs  to  the  system 
— o 

A,  B  B  ,  C  ,  C  ,  C  =  constant  matrices  of  respective  suitable 

—  — u.  — v.  — v  - ■ w  — o 

l  l 

dimensions 

If  we  want  the  plant  to  operate  at  state  vector  and  suppose  that  this 
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state  can  be  made  an  equilibrium  state  of  the  system,  i.e, 


N 


N 


x  =  [I  -  A  - 


u 

E 


v 


u  .  ■,  B 

.  n  id  — u . 
1=1  1 


,-1 


E  Vid  -V1  {  [C  C  ]x 

1  =  1  1  —V  — W 


^d 


^d 


+  C  } 
— o 


(V . 12) 


The  above  expression  is  obtained  by  solving  equation  (V.ll)  with  the 

left  hand  side  set  to  zero  and  x_(k)  ,  u.(k)  ,  w(k)  ,  _v(k)  respectively 

replaced  by  x_,  u,,  w,,  v.  u,,  w,,  v,  are  control  inputs  to  keep  the 

— d  — d  — d  — d  *  —  d  — d  — d 

system  at  equilibrium. 


Define  the  error  state  variables  as  follows: 

6X(k)  =  X(k)  -  XJ 
—  —  — d 

or  X (K)  =  Xd  +  6X(k) 

Substituting  (V.13)  into  equation  (V.ll),  we  obtain: 


(V . 13) 


N 


N 


u 


V 


fix (k+1)  =  [A  +  Z  u.(k)  B  +  Z  v.(k)  B  ]  6X(k)  + 
—  —  l  — u.  .  ,  i  -at.  — 

i=l  i  i=l  i 


N 


N 


u 


V 


[  Z  (u .  (k) 
i=l  1 


u.  J  B  +  Z  (v .  (k) 
id  — u.  .  ,  l 

l  1=1 


v.  J  B  ]  X  + 
l  d  — ^ v  .  ~  d 

l 


[C  C  ] 
—v  w 


X  -  Y.d 


w  -  w. 
—  — d 


(V.14) 


N 


N 


u 


v 


Let  A  =  A  +  Z  u  B  +  Z  v  B 

— d  —  .  ,  id  — u.  id  v . 

1=1  l  1=1  l 
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6u  .  (k)  =  u  .  (k)  -  u  .  , 

1  1  id 

6v .  (k)  =  v .  (k)  -  v . 
i  l  id 

6w .  (k)  =  w .  (k)  -  w 
l  l  id 

Equation  (V.14)  can  be  rewritten  as: 


N 


u 


N 


y 


6X(k+l)  =  [A  +  £  6u.(k)  B  +  £  6v.(k)  B  ]  6X(k)  + 

—  ~  d  .  n  l  — u .  .  l  — 

i=l  i  i=l 


D 


5u  (k) 
5v(k) 
6w(k) 


(V . 15) 


where. 


D  = 

[D  | 

D 

|d  ] 

u 

V 

1  w 

D  = 

[B  • 

B  x  I  .  .  . 

.  1  B  x 

— u 

— d 

1  u  — d 

U 

D  = 

[B  . 

X  j 

+  C  .  . 

B  .  x,  • 

— V 

~V1 

d 

^1 

“d 

V 

N 


v 


D  =  C 

— W  “V 


Hence,  a  chosen  /  0  generally  introduced  additive  control-variation 
terms  for  all  the  multiplicative  control  variables,  and  there  is  no 
"separation"  of  additive/multiplicative  control  inputs  (except  that 
w,  if  present,  remains  only  additive). 

If  _dX(o)  is  the  initial  error  state,  and  if  we  can  find  a  control 
sequence  of  length  N  such  that  6X(N)  =  o_,  then  defining 


M. 
“ J 


[Ad  + 


N 

u 

£ 

i=l 


6u.  (j)  B 

l  — 


N 

v 

+  £  6v  (j)  B  ] 

.  i  -v. 

i=l  l 


l 


(V. 16) 
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From  (V.15),  we  can  write 


6X(1)  =  M^6_X  (o)  +  D 


6u  (o) 
6v(o) 
6w  (o) 


6X(2)  =  M^XU)  +  D 


6u(l) 
6v(l) 
5w  ( 1 ) 


(V . 17a) 


_6X(N)  =  0 


<5X(N-1)  +  D 


6u(N-l) 

6v(N-l) 

5w(N-l) 


or,  eliminating  all  _6X(j)  for  j  >  o,  we  obtain: 

-IVi  V2  •  •  •  V  ^(o)  “  ID  I  «N-1  £  I  •  •  -  I  ®N-1  •••«!>  £1  x 


[^(N-i)  6vT(N-1)  6wT(N-1)]T 
[6uT(o)  6vT(o)  6wT(o)]T 


(V. 17b) 


Thus  the  problem  is  to 

a)  find  a  sequence  of  controls  which  satisfy  the  above  equation, 
and 

b)  optimize,  in  some  sinese,  the  performance  of  the  system. 
Condition  (a)  is  the  controllability  problem,  and  we  shall  show  that 
this  has  solution,  in  general;  even  if  only  one  control  variable  is 
available  and  it  is  both  additive  and  multiplicative.  The  more  difficult 
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question  is:  how  to  use  the  freedom  of  the  richer  control  structure 
to  optimize  performance  (b)? 

Suppose  that  we  would  like  to  minimize  the  vector  length  of  the 
solution  control  sequence  appearing  on  the  right  hand  side  of  the 
equation  (V.17),  i.e.. 


N-l 

ii  ,  1 2  ii  ii2  ii  .  |  ?  1/9 

min  [  E  (a  |  |_6u(i)  |  |  +  a  |  |  5v(i)  |  |  +  a  |  6w(i)  |  )  ] 

.  ,  u  V  w 

1=1 


Assuming  that,  minimization  of  the  above  expression  is  directly  related 
to  the  elapsed  time  for  the  control  efforts  to  bring  the  plant  back  to 
its  equilibrium  state,  in  other  words,  the  faster  the  system  reaches  its 
equilibrium,  the  less  power  it  consumes  (in  a  physical  system,  _6u,  5v, 

6w  represents  power  consumed) . 

Pick  one  of  the  control  variables,  say  6v^,  and  initially  set  it  to 

zero 


6v1(j)  =  0 


j  =0,  1 ,  . . . ,  N-l 


Then  use  the  remaining  control  variables  to  establish  a  set  of  state- 

transition  matrices  M  ,  . ..,  M  such  that  each  M.  has  a  "fast"  eigen- 

o  N-l  3 

vector  m.  (i.e.,  M.  m.  =  z.  m.  with  I  z .  <<  1 )  and  the  set  (m  ,  in  ,  ..., 
-J  “J  J  1  “J  J  -o  1 

m^  are  as  nearly  as  orthogonal  as  possible  (assuming  that  N  =  n  = 
dimension  of  the  state  space)  i.e.,  we  want  them  to  be  a  basis  in  the 


state  space  _T 

N  N 

u  v 

where  M.  =  [A  +  E  6u.(j)  B  +  E  6v.(j)  B  ] 

— j  d  .  l  — u .  .  „  1  —17  , 


(V.18) 


i=l 


i=2 


The  idea  behind  the  above  algorithm  is  that,  to  minimize  the  control- 
variation  efforts,  we  would  like  to  have  6v  as  small  as  possible,  in 
some  sense.  This  is  achieved  if  such  a  set  of  "fast"  eigenvectors  can 


be  synthesized  using  the  other  control  variables  (in  particular,  the 
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the  purely  multiplicative  ones,  the  u's,  may  not  involve  much  energy 
comsumption  so  their  magnitudes  need  not  be  minimized.  This  is  true 
for  a  solar-assisted  heat  pump  system).  Thus,  if  6v^ can  be  kept  small, 
the  perturbations  caused  by  6v^  0  on  the  eigenvectors  m.  can  be  small 
so  the  synthesis  will  still  be  fast  to  a  close  approximation. 

To  implement  the  algorithim  as  closed-loop  control  switching, 
the  current  measured  state  vector  would  be  analysized  in  the  basis  of 
the  fast  eigenvectors,  e.g.: 


<5X 
— o 


y  m  +  y  m  +  ...  +  y  m 
o — o  1—1  n- 1 — n- 1 


(V. 19) 


Identify  the  largest  component,  let  us  say 

y.  >  y  all  j  4-  i 
1  -  1 

Adjust  the  "other"  control  variables  to  produce  state  transition  matrix 

M  . 

— i 

Then,  after  one  sample  period,  if  no  further  state  disturbance 
occurs,  the  next  state  would  be: 


6X  =  Mi  6Xq  +  D 


6u  (o) 
5v  (o) 
6w  (o ) 


n-1  « 

Z  y  M  m 
j-o  3  1  J 


+  y  .  z  .m  .  +  D 

l  i—i  — u 


6u  +  D'  6v’  +  D  6w  + 
—  —  v  —  ~ w 


(B 


x  i 

V1  ”d 


+  C  )  6 


v. 


V, 


D’  =  {last  (N  -1)  columns  of  D  } 
— v  v  —v 


(V . 20) 


where 
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6v *  =  {last  (N  -1)  components  of  6v} 

Since  |z_^|<<1,  the  dominant  np  component  of  6X  has  been  greatly  reduced 
by  the  action  along  the  "fast"  eigenvector  direction. 

The  additive  terms  in  _6u_,  _6w,  and  6v^  (j  ^  1)  may  be  viewed  as 
constant  forcing  inputs  using  the  sample  period,  causing  some  further 
incremental  displacement  of  the  state.  The  remaining  additive  control 
term  6v^  should  be  set  at  the  appropriate  value  for  "deadbeat  control". 
For  example,  suppose  we  ignore  the  other  additive  terms  for  the  present, 
and  set 


6X,  =  M.  SX  +  V<5v, 

— 1  —i  — o  —  1 


Then  SX  ->  o_  in  one  sample  period  if  the  equation  o^  =  M  6Xq  +  _6v^  has 

~-l  r 

a  solution  6v^(o).  This  is  possible  if  and  only  if  5X^  =  FL  _  for 

—  1  V 

some  scales  a  i.e.,  _SX^  happens  to  lie  along  the  vector  _  in  state 
space.  If  this  is  the  case,  then  choose  6v^  equal  to  -a. 

Likewise,  6X^  ■+  _o  in  two  sample  periods  with  an  open-loop  control 
sequence  {6v^(o),  6v^(l)}  if  and  only  if  the  following  equation  has  a 
solution  sequence 


£  =  Mi(1)  6X(1)  +  l6v1(l) 


or  c>  =  6X(o)  +  ^  6v1(o)]  +  ^  Sv^l) 


The  above  equation  has  solution  if  and  only  if 


6X(o)  =  otM*  .  V  ,  m.)  s  M.^  V 
—  X ( o )  -  +  — i(o)  — l(l)  — 


for  some  asj6  the  two-step  open  loop  sequence  would  be  Sv^(o)  -  -a 

Sv^l)  =  -3. 
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In  general,  for  an  n-step  dead-beat  control  sequence,  _6X(o)  has 
to  be  in  the  following  direction: 


n-1  1 

6X(o)  =  oM.  (  ,  v  +  E  [  tt  M . X.  V 

~  "l(o)  “  1=1  k=o  l(k)  1  “ 


(V . 21) 


where  6v^(o)  =  -a 

and  Sv^(l)  = 


1  =  1 ,  . . . ,  n-1 


Note,  however,  that  only  6v^(o)  would  be  used  in  closed-loop 
control  (eqn.  (V.20)),  after  one  sample  period  a  new  feed-back  state- 
measurement  would  be  available.  Nonetheless,  all  B’s  have  to  be 
calculated  in  order  to  find  a. 

The  other  additive  forcing  terms  do  not  really  complicate  the 
matter  very  much.  If  we  have,  for  example, 


o  =  M.,  .  5X(o)  +  ^6v_ (o)  +  ^  6w(o) 
—  —  i  ( o )  —  —  1  — w 


(V.22) 


we  must  have 


6X(o)  =  aM*  v  v  -  M*  v  6w(o) 
—  1  ( O )  —  —1  ( o )  — w 


then  6v^(o)  =  -a  is  chosen,  here  we  have  the  freedom  of  adjusting 

6w(o)  for  a  proper  a.  However,  if  6u(o)  takes  the  place  of  6w(o) 

in  (V.22),  a  cannot  be  varied  because  6u(o)  has  been  fixed  by  the 

process  of  establishing  M.  ,  . 

— i(o)  . 

From  a  slightly  different  point  of  view,  equation  (V.17)  defines 
an  open-loop  sequence  of  control  such  that: 


-[M.  ,  ..  M  , 

— i(n-l)  — i(n-2) 


M.(o)]«X(o)  -  [D  |  D 
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M.  .  ...  M  D]  x 

i (n-1)  -i(l)  — 


[ 6uT(n-l)  6vT(n-l)  6wT(n-l)]T 


[6uT(o)  6vT(o)  6wT(o)]T 


(V . 23a) 


Here  M^_^,  ^N-2’  etc*  -*-n  eTn'  (v*17)  are  replaced  by  M-j_  (n_j)  » 

M.  ,  .  ,  etc.  . . . 

—i  (n-2) 

Assuming  that  6v  (j)  (j  =1,  ...,  n)  is  negligible. 

Equation  (V.23a)  can  be  rewritten  as  (compare (V . 16)  and  (V.18)): 


-[M  M  ...  M  (o)]5X(o)  =  [D  |  M  D|  ... 

— i(n-l)  — i(n-z)  — i  —  —  1  — i(n-l)— 1 


M.  ,  .  ...  M.  ,  .  D]  x 

l (n-1)  — i(l)  - 


[6u^(n-l)  6v^(n-l)  6w^(n-l)]^ 


T  T  T  T 

[6u  (o)  6v  (o)  6w  (o)] 


(V . 23b) 


i(o),  i(l),  etc.  ...  denote  the  index  of  the  particular  "fast  eigenvector 

matrix"  that  would  be  selected  at  each  sampling  period.  That  is,  in 

the  first  sample,  solving  equation  (V.23)  give  us  6v^(o)  (we  obtain 

6v^(j),  j  =  0,  ...,  n-1,  but  only  6v^(o)  is  needed  in  this  first  step) 

with  all  other  controls  having  been  given  in  the  process  of  establishing 

M,  , N,  M  ,  „N ,  .  .  .  M  .  .  .  Using  equation  (V.17a),  the  state  vector 
— o(n-l)  — o(n-2)  — o(o) 

at  the  second  step  is  calculated: 


6X(1) 


M  SX(o)  +  D 
— o  —  — 


6u  (o) 
6v  (o) 
6w(o) 


where  M  =  M  .  .  ,  B  6v, (o) 

— o  — o(o)  +  v^  1 


(V.17a) 
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The  control  sequence  length  now  is  (n-1) ,  solve  for  6v  (1)  by  using 
the  modified  equation  (V.23b) 


[-l(n-l)  — 1  (n-2)  -1(1)]— (1)  ^  I  — i (n-l)  -  I  ** 


Ml(n-1)  *  *  *  -1(2)"  X 


[6u^(n-l)  Sv^(n-l)  6w^(n-l) ] ^ 


T  T  T  T 

[  6u  (1)  6v  (1)  6w  (1)  ] 


The  state  vector  at  the  third  sample  is 


6X(2)  =  M  6X(1)  +  D 


6w(l) 

<Sv  (1) 

5w(l) 


where 


=  M 


KD 


B 

— v. 


(1) 


And  then  repeat  the  process  by  using  equation  (V.23b)  to  calculate  the 
control  6v^(2),  etc.  ...  From  the  equation  (V.23b),  we  can  write: 


-_6X(o)  =  £ 


<5v^  (n-l) 

6v  (n-2) 

• 

6v1(o) 


+  (a  vector  independent  of  6v^)  (V.24) 


where  d>  =  [M  ^  ...  M  ^  .  D  |...|m.^  n  D  ] 

-  -i(o)  — i(l)  — i(n-l)  — v  1  i ( o )  ~v1 

D  =  the  part  of  D  that  multiplies  dv.,  (it  was  called 
r  — V  1 

in  equation  V.21). 


D  =  B  X,  +  C 

~vi  _d  “'i 


V 
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The  question  is  what  properties  the  matrix  £  should  have  in  order 


to  minimize  the  solution  [6v^(n-l)  6v^(n-2) 


6v  (o)]  ?  Since  we 


are  inverting  this  matrix  to  solve  for  the  6v  ’  s ,  the  matrix  should  be 
nonsingular  and  well-conditioned.  That  is,  the  columns  of  cji  must  be 
linearly  independent,  should  be  as  nearly  orthogonal  as  possible,  and 
£  should  be  uniformly  large  in  some  sense  in  order  for  the  five's  to  be 
minimized.  In  the  sense  of  the  spectral  norm,  the  eigenvalues  of  tj> 
should  all  be  large  in  magnitude. 

In  general,  these  conditions  may  perhaps  all  be  met,  as  will  be 
proven  in  the  following,  if  and  only  if  ,  etc.  each  has  one 

"fast"  eigenvector  and  these  eigenvectors  are  nearly  orthogonal. 

The  sufficient  condition  of  the  above  statement  will  be  proven 
as  follows,  for  the  necessary  condition  it  is  similar. 

Let  the  last  column  of  6  be  i  ,  it  can  be  written  as: 

-  — n 


<b  =  M .  ^  x  D  or  D  =  M .  ,  .  cb 
—  -1(0)  —v1  v  —l(o)  -^n 


b  can  be  resolved  in  the  basis  of  the  fast  eigenvectors  {m.  ,  , 
■ti  —  x(o). 


m 


-i(l) 


,  m 


i (n-1) 


} 


or 


(p  =  a  .  m .  ,  .  +  a  m.  a  .m.,  ,  N 

— n  n,0  —1(0)  n,l  — i(l)  n,n-l  — i(n-l) 


M.M  i  =  a  M  m  +a  M  m  + 
— x(o)  n,0  —1(0)  —1(0)  n,l  —1(0)  — i(l) 


D  =  a  nz .  ,  .  m  +  a  ,  M.,  >  1.,  ,  +  ... 

— v^  n,0  1(0)  —1(0)  n,l  —1(0)  — i(l) 


Because  D  is  a  constant  vector,  and  z.,  x  is  the  fast  eigenvalue,  hence 
"  1(0) 


an , 0  i(o) 


constant 
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is  very  small,  thus  ^  must  be  large.  Therefore,  the  component 
of  along  the  direction  m.  ^  ^  must  be  large. 

Likewise,  the  next-to-last  column  of  is: 


<}>  ,  =  m.-J:  .  mT^x  d 

— n-l  — 1(0)  — i(l)  — v. 


or , 


M  .M  (f>  =  D 

— i(l)  — 1(0)  —  n-l  — v. 


Similarly  to  the  previous  case,  the  component  of  vector  M. ,  .  d> 

r  — 1(0)  -"n-l 

along  the  direction  m.  ,  N  must  be  large,  or  4'  must  have  a  large 

— 1  (,  1 ;  — : n—  1 

component  in  the  direction  of  vector  M.^  .  m.,,.  which  is  not  close  to 

-i(o)  — i(l) 

the  direction  of  m.,  N.  Suppose  it  were,  we  can  write: 

— i(o) 


~-l 

M  m.  . 

— i(o)  — i(l) 


am 


i  (o) 


or 


m  -  aM  m  =  az .  .  .  m 
-i(l)  — i(o)  — (o)  i(o)  — x(o) 


So,  vectors  m.  ,,  .  and  m.  .  .  are  of  close  direction,  this  contradicts 
-i(l)  — i(o) 

the  assumption  that  all  the  eigenvectors  m.  . ..,  m.  .  .  ,  (j,  k  =  0,  ...,  n-l) 

x  J  )  i  (k ) 


are  nearly  orthogonal, 


-1 


Therefore,  direction  M.  .  .  m..,.  is  not  close  to  direction  m.,  ,  , 

— x(o)  —  i(l)  — i(o) 

matrix  ^  can  be  thought  of  as  amplifying  the  component  of  m^^ 


that  is  in  the  direction  of  m^^,  but  that  component  is  very  small 

because  m.,,x  and  m. ,  s  are  nearly-orthogonal .  In  other  words,  column 
-i(l)  i ( o ) 


(j>  is  nearly-orthogonal  to  col 
column  of  <j)  is: 


umn 


y  -t 

n-l 


Continuing,  the  second-to-last 


~-l  ~-l  ~-l 

6  =  M.  M.  ,  .  M.  .  D 

— n-2  — i(o)  — i(l)  — 1(2)  v . 
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or  M  M  M  s 

-i(2)  -i(l)  i (o ) 


4-2 


=  D 
— v. 


Hence,  Mq  ,  ^  -^-2  ^aS  a  ^ar§e  component  along  the  direction  of 

m. .  .  or  4  _  has  a  large  component  along  the  direction  of  M,,  N  M.^.N 

— 1(2)  -nv-2  — i(o)  —i(l) 

m./ON.  M. -  N  and  respectively  amplify  the  components  of  m. 

-i(2)  i ( o )  — i(l)  — i (2) 

that  are  in  the  direction  of  m.,  ,  and  m.,,N  respectively,  but  these 

i (o )  — i(l) 

components  are  small  because  m. is  nearly  orthogonal  to  m.  ,  ,  and 

-i(2)  J  -i(o) 

m..,x,  therefore,  direction  M.^  N  M.'L.  m.,„.  is  not  close  to  the 
-i(l)  — i(o)  -i(l)  — i (2) 

direction  of  in.,  *  and  m.,^N,  or  the  direction  of  the  second- to-last 
i ( o )  -i(l) 

column  of  h  i  is  not  close  to  the  directions  of  <j>  .  and  <j)  , 

-  — n-2  — n-1  — n 

and  is  nearly-orthogonal  to  ft  and  cp  (because  components  of  6 

— n-1  — n  n 

along  the  direction  of  (p  and  6  ^  are  very  small)  . 

Similarly,  we  can  prove  that  all  the  directions  of  all  the  columns 
of  <£>  are  not  close  to  one  another,  or  they  are  linearly  independent, 
and  are  relatively  orthogonal  to  one  another. 

Let  U  be  a  set  of  basis  rectors  which  are  in  the  direction  of 
m.,  M .  \  .  m.  ,  .  ,  M .  ^  M.  ,  .  m  ,  etc.  ...,  respectively. 

Because  these  vectors  are  nearly  orthogonal,  if  we  express  columns 
of  <£  in  this  basis,  £  will  be  composed  of  mainly  diagonal  elements. 
These  diagonal  elements  of  <£  will  be  its  component  in  the  direction  of 
the  basis  vectors  U  and  these  elements  are  large,  hence,  all  of  the 
eigenvalues  of  £  are  large  or  £  is  uniformly  large  in  terms  of  spectral 


norm. 


Solving  equation  (  V.24),  we  obtain  the  entire  deadbeat  open-loop 

sequence.  But  the  practical  problem  of  finding  the  control  sequence 
6v_^  remains  difficult,  since,  as  we  see  previously,  the  matrix  depends 
on  the  matrices  (i  =  1,  n;  j  =  1,  ...,  n)  ,  which,  in  turn, 


I-6HI 
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changes  at  every  sampling  step.  So  this  would  mean  inverting  a  different 
at  each  sample  and  this  is  computationally  infeasible  with  micro¬ 
processor  . 

Until  now  we  have  assumed  that  null-controllability  of  the  system 
is  possible  with  the  two-step  process  of: 

a)  First  set  5v  (j)  =0  j  =  0,  1,  n;  then  find  n  different 

sets  of  values  for  all  the  other  fiu's,  fiv’s  that  give  n  relatively 
orthogonal  fast  eigenvectors  for  the  corresponding  state-transition 
matrices . 

b)  Then  find  6v^  values  that  tend  toward  deadbeat  control  closed 
loop,  assuming  that  these  values  are  small  and  do  not  significantly 
affect  the  fast  eigenvectors  found  in  (a) . 

In  the  following,  we  will  examine  this  assumption  to  see  whether  it  is 
indeed  possible  to  have  null-controllability  when  fiv^  in  fact  appears 
both  multiplicatively  and  additively. 

First  consider  the  case  n  -  2,  assuming  that  _5X(2)  =  _0  can  be 
achieved.  From  eqn.  (V.17a),  we  can  write: 

6X (1)  =  M.  .  . 5X(o)  +  D 
—  — x  (o) —  - 

and 

5X (2)  =  0  =  Mi(1)5X(l)  +  D 


6u(l) 

6v(l) 

fiw(l) 


fiu  (o) 
fiv  (o) 
6w(o) 


or 

6w(o) 

fiu(l) 

0  =  M.  ...  [M.,  .  5X (o)  +  D 
-  -i(l)  — i(o) — 

fiv  (o) 

]  +  D 

fiv(l) 

fiw(o) 

5w  (1 ) 

_  _ 
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or 


0  =  M .  .  . M .  ,  5X(o)  +  M .  .  . D  6vn  (o)  +  D  dvn (1)  + 

— i(l)— i(o) —  — i(l)— 1  — v  1 


(a  vector  independent  of  6v^(o)) 


(V . 25) 


compare  equations  (V.16)  and  (V.18),  we  obtain 


M.  =  M-  +  B  6vn  <k) 
— i(k)  -i(k)  v  1 


(V. 26) 


Substituting  (V.26)  into  (V.25) 


0  =  [M.(1)  +Bv^Vl(l)]  [M.(o)  +  Bv^6v1(o)]6X(o)  + 

[M  +  B  <5v  (1)]  D  dv  (o)  +  D  dv  (1)  + 

— i(l)  -^v^  1  — 1  — 1 


(a  vector  independent  of  dv  (o)). 


or 


0  =  -dX(o)  -  [M.  ,  s  +  B  dv  (o)]  ^  D  dv  (o) 
—  -  —  i(o)  — v^  1  — v  1 


-[M.  +  B  dv  (o)  ]  1  [M  ,  ,  +  B  5w  (1)]  XD  dv.  (1)  + 

— l(o)  — V^  1  —  l(l)  — v^  1  v^  1 


(a  vector  dependent  of  du’s,  dv’s,  dw's  and  other  than 


dv1(o),  dv1(l)) 


(V . 27) 


(V.27)  is  obtained  by  assuming  that  dv^(o)  and  dv  (1)  are  very  small  and 
can  be  neglected  in  the  last  term  of  the  equation. 

Equation  (V.27)  is  true  if  and  only  if  we  can  find  dv^^  and 
dv  (1)  such  that  the  vector  function 


C  =  [M .  ,  x  +  B  dv  (o)  ]  dv  (o)  +  [M.  ,  .  +  B  dv  (o)  ] 

—  —  i(o)  V-^  1  — v^  1  —  i(o)  — v^  1 


x 
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[Mi(l)  +BVi6v1(1)]-1Dv^v1(1) 


(V . 28) 


covers  the  entire  plane,  i.e.,  values  of  6v^(o)  and  6v  (1)  can  always 
be  found  such  that  this  function  equals  the  arbitrary  vector  given  by 
the  other  terms  in  equation  (V.27).  C  can  be  written  as  the  sum  of  the 
two  vector  functions: 


C  =  sj  Sv^o)  ]  +  t_[6v1(o),  6v1(l)] 


(V . 29) 


The  locus  of  _s  will  be  examined  first.  _s  is  a  curve  in  R^,  passing 

through  the  origin  when  6v^(o)  =  0  and  asymptotic  to  the  straight  line 
~  4 

Hv  6v^(o)  as  <5v^  ( ° )  -*•  0.  For  6v^(o)  ^  0,  there  may  be  up  to  two 

finite,  real  values  of  6v_,  (o)  for  which  M.  .  .  +  B  .  6v_  (o)  is  singular. 

1  — i(o)  — v^  1 

Let  v  *  be  such  a  value.  Then  as  6v_  (o)  ->  v*,  s[6vn(o)]  °°  asymptotic 

o  1  o  —  1 

to  the  null  space  of  [M. ,  N  +  B  v  *],  which  is  denoted  as  N  *  (Fig.  V.l). 

— i(o)  — v^  1  o 

If  there  are  two  singular  values  v  *  and  v  **  with,  say  0  <  v  *  <  v  ** 

o  1  o  o 

then  the  locus  consists  of  two  curves,  each  tending  to  00  as  6v^(o)  v^* 


and  v  **,  respectively  (Fig.  V.2). 


o 


-1 


Assuming  B  is  nonsingular,  as  6v.  (o)  ±°°,  s  ->  B  D  ,  if 

vi  1  "  ~^i  ^1 

B  is  singular: 

^1 


s  =  [M. ,  .  +  B  6v  (o) ]  ^  D  6v_ (o) 
-1(0)  — v  1  -v  1 


(V . 30) 


or  [M.  ,  x  +  B  <5v-(o)]s  =  D  6v  (o) 
—1(0)  — v  1  —  -^1  1 


when  6v^(o)  -*  °°,  s  ->  °°,  asymptotic  to  the  null  space  of  B^  .  Does  the 
curve  s  loop  and  cut  itself?  I.e.,  can  we  have  some  _s  for  which 


s=[M./n+B  a]^D  (a)=[M.,x+B  8]  D  8 

-1(0)  — v  — -1(0)  vx  -v1 


(V.31) 
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for  a  /  g,  a  <  v  *,  3  >  v  *,  a  and  6  are  two  distinct  values  of  6v_,  (o)  . 

o  o  1 

(v.31)  can  be  rewritten  as: 


[M.  .  .  +  B  a]s  -  D  a  =  0  =  [M..  +  B  3]s  -  D  6 

-i(o)  -v1  -  -v1  i(o)  — -  — v 


or  (a  -  3)  (B  s  -  D  )  =  0 
— v  —  — V  — 

1  1 


Since  a^3 


B  s  -  D  =0 

T  —  — V  — 

1  1 


This  implies 


M.  ,  .  s  =  0  (V . 32) 

— i(o)— 

But  this  is  impossible,  since  M.  ,  N  is  a  state-transition  matrix  for 

— i(o) 

a  discrete-time  system  and  it  cannot  be  singular.  The  only  remaining 

possibility  then  is  _s  =  _0  as  a  point  of  intersection,  but  from  (V.31) 

this  leads  to  [M..  .  +  B  a]  ^  and  [M..  N  +  B  3]  ^  are  singular.  This 

i(o)  — v  — i(o)  — v 

is  a  contradiction  since  [M.  .  x  4-  B  a]  and  [M.  ,  N  +  B  3]  are  both 

-i(o)  — v  ~i(o)  v x 

well-defined,  finite  matrices  which  are  nonsingular  and  have  nonsingular 
inverses  [M.  .  N  +  B  a]  ^  and  [M. ,  ,  +  B  3]  ^ .  Therefore,  s  =  0  is 
not  a  point  of  intersection,  either.  In  short,  the  curve  does  not  loop 
and  cut  itself. 

Now  consider  the  function  t 


t[6v  (o) ,  6v  (1)  ]  =  [M.  (  ,  +  B  6V  (o)  ]  1  [M  ,  +  B  6v  (1)  ] 

—  1  1  — l(o)  — V ^  1  —1(1)  1 


D  6v  (1) 
1 


(V .  33) 


for  each  fixed  value  of  6v^(o),  _t  is  a  function  of  5v^(l).  t_  can  be 
thought  of  as  a  vector  displacement  which  is  added  to  sJSv^(o)]  for  each 
fixed  value  of  6v^(o),  _t  sweeps  out  a  curve  in  the  plane  which  passes 
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through  the  point  sJSv^Co)]  for  6v  (1)  =  0.  This  curve  goes  to  °°  as 


(1)  ->  v  *  or  v  **,  where  det[M.  ,  .  +  B  6v  (1)]  =  0,  the  curve  goes 
1  1  i  — i(l)  — v^  1 


to  infinity  asymptotic  to  the  null  spaces  N  *,  N  **  of  fM  ,  „  +  B  v  *lx 

—1—1  — i(l)  -v  1 

[-i(o)  Svl(0))  or  [^i(l)  +^1V1**1  [“i(o)  +£VlSvl(o)1- 


,-1 


If  B  is  singular,  t  00  asymptotic  to  the  line  til.  ,  .  +  B  v„  (o)  ]  \ 

— —  —1(0)  — v  1 


[null  space  of  B_  ]  when  6v  (1)  ->  ±°°. 

V1  1 


•1 


If  B  is  nonsingular,  when  6v,  (1)  ±°°,  t  [M.  .  .  +  B  v,  (o )  ]  x 

— v  1  —  —1(0)  v  1 

-1 

B  D  . 

— V  — V 
1  1 

From  an  exampled  figure  in  Fig.  V.3,  we  see  that  each  branch  of 
sj6v^(o)]  has,  attached  to  each  point  on  it,  a  displacement  curve 
_t[6v^(o),  6v^(l)]  with  the  corresponding  value  of  6v  (o)  ,  sweeping  off 
to  00  as  a  function  of  6v^(l). 

It  seems  quite  likely  that  these  curves  will  cover  the  whole  space 

2 

|R  ,  in  other  words,  there  always  exists  a  null-control  sequence  (6v^(o), 

6v  (1)}  for  an  arbitrary  "error"  state  vector.  These  curves  certainly 

2 

cover  a  neighborhood  of  the  origin  in  IR  ,  since  they  are  asymptotic  to 
straight  lines  there.  Here  the  relatively  small  values  of  6v  (o) , 

6v^(l)  imply  only  small  curvature  will  be  introduced  by  the  6v^  ^  0 
effect,  as  a  multiplicative  control  variable,  on  the  state-transition 
matrix . 

But  this  picture  also  raises  the  question  of  how  one  could  compute 


the  accurate  values  of  6v^  for  open-loop  deadbeat  control.  The  previous 

algorithm  totally  ignores  the  curvature  of  the  covering  loci  and  would 

2 

be  accurate  only  in  the  limit  as  6v^  -*  0  near  the  origin  in  IR  (in  eqn. 

(  V.18)  if  6v  (o)  and  B  6v  (1)  are  ignored,  s_  and  _t  will  be  straight 

V1  ~^1 

lines  and  are  linear  functions  of  6v^(o)  and  6v^(l)  respectively, 
that  is  exactly  how  state- transition  matrix  M.  is  set  up  in  equation 
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(  V.18)).  That  algorithm  also  involved  computation  of  £  b6X(o)  (or 

at  least  the  bottom  rows  of  all  possible  matrices)  either  on-line  or 

off-line  (with  bottom-row  storage),  which  could  be  quite  memory-demanding. 

An  alternative  that  would  account  for  the  curvature  effects  as  well, 

would  be  to  store  a  grid-map  of  pairs  of  values  [6v^(o),  6v  (1)]  that 

2 

cover  a  sufficiently  large  neighborhood  of  the  origin  in  |R  ,  then  simply 

use  the  current  measured  value  of  _6X(o)  to  look  up  the  correct  current 

control  value.  Actually,  only  6v^(o)  values  are  needed  with  closed 

loop  control,  but  they  are  needed  for  a  two-dimensional  (n-dimensional ,  in 

2 

general)  grid  of  points  in  |R  .  This  can  become  equally  demanding  of 
memory.  However,  since  if  the  grid  has  jd  points  along  each  dimension 
inlRnspace,  the  grid  has  j}n  points,  hence,  we  need  p>n  words  of  memory 
to  store  6v^(o)  values.  The  previous  algorithm  requires  n  x  n!  words 
to  store  all  possible  bottom  rows  of  cji .  Using  Stirling's  formula: 

„  (n+1)  -n  fz —  m>n  jz — 

n.n!  =  n  e  v27Tn  =  n(— )  /2Trn 

e 

we  see  that  this  will  grow  even  faster  than  j)n  for  n  large. 

For  n  >  2,  the  loci  covering  arguments  presented  above  can  be 
extended  in  an  obvious  way  to  cover  a  neighborhood  of  the  origin  in  Rn 
with  "swept  out"  manifolds  of  growing  dimension  by  families  of  parametric 
curves 

s[  ^  (o)  ] 

sa  +  t[6via(o),  6v^ (1) ] 

_sa  +  +  r[6via(o),  6v^b  (1) ,  6v1(2)] 


etc . 
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Where  s_a  is  a  fixed  point  on  the  s_(6v  (o)  ]  curve,  s3  +  jtb  is  a  fixed 
point  on  the  s_  +  _t  surface  and  the  family  of  all  _r  displacements 
passing  through  all  such  fixed  points  covers  a  three-dimensional 
manifold,  etc. 

It  is  clear  that  these  manifolds  will  always  cover  a  neighborhood 
of  the  origin  in  Rn,  since  all  the  curves  are  asymptotic  to  the  linear 
subspaces : 


M.)  D 
-1(0)  v 

MT)  ,  M-l  D  +  D 

— l(o)  —1(1)  —v^  — l(o)  — V 


'"-1  ~-l  ~-l  ~-l  ~-l  ~-l 

M .  \  .  M  M  D  +  M.  ,  M . .  D  +  M.,  .  D 

— i(o)  — i(l)  — x (2)  -v  — i(o)  -i(l)  — i(o)  —v 


etc . 

This  is  assured  by  the  fact  that  the  M  matrices  are  designed  to  have 

relatively-orthogonal  fast  eigenvectors  which,  for  example,  make 
~-l  ~-l  ~-l 

M.  .  N  D  and  M.,  x  M./1N  D  not  parallel  or  nearly  parallel. 

— i(o)  — v  — i(o)  — i(l)  — v 

Basically,  the  deadbeat  control  approach,  that  has  been  talked 
about  until  now,  is  that  assuming  one  of  the  multiplicative  additive 
control  inputs,  say  6v^,  is  negligibly  small  so  that  a  set  of  fast 
eigenvectors  can  be  approximately  obtained,  and  these  fast  eigenvectors 
are  used  to  restore  the  system  state  back  to  its  equilibrium  point. 

The  solutions  to  the  problem  of  finding  6v^  are  proposed,  five’s  are 
obtained  either  by  inverting  a  matrix  t|>,  or  by  setting  up  a  map  of 
values  of  five's  in  a  multi-dimensional  space. 

Although  this  method  looks  promising,  there  are  several  problems, 
as  follows,  that  remain  to  be  solved.  The  global  controllability  has 
not  been  examined  and  it  is  doubtful  that  based  on  the  richness  of  the 


system  structure,  we  can  obtain  global  controllability  for  a  general 


_ 
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bilinear  system.  In  appendix  IV,  a  particular  bilinear  system  (the 
solar-assisted  heat  pump  system)  is  examined  in  terms  of  its  eigen¬ 
values.  As  we  see,  the  eigenvalues  of  this  system  are  always  negative, 
and  this  violates  the  sufficient  conditions  for  complete  controllability 
[2],  although  it  is  not  enough  for  us  to  draw  any  conclusion  yet,  but 
it  appears  the  system  is  hardly  globally  controllable. 

The  property  of  the  set  of  the  equilibrium  points  (eqn.  V.12) 
is  not  investigated  either.  This  should  be  done  and  hopefully,  to  be 
worthwhile,  the  equilibrium  set  will  cover  the  whole  state-space  or 
at  least  a  large  portion  of  it. 

There  also  exists  the  practical  problem  of  implementing  the 
controller.  This  approach  can  be  referred  to  as  supoptimal,  never¬ 
theless,  nothing  has  been  done  to  prove  that  it  is  nearly  optimal. 
Neither  is  there  any  proof  that  we  can  always  obtain  a  set  of  fast 
eigenvectors . 

Due  to  the  time  limitation  we  have  not  been  able  to  address  every 
aspect  of  this  approach,  however,  we  believe  that  this  approach  could 
be  feasible  especially  in  solar-assisted  systems. 
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Fig.  V.1:  Locus  of  s  (one  value 
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Locus  of  s. 
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[Mkoj+B*.  Sv^(0)]  Mi(,)DVi 

Pig.  v . 4 :  Locus  of  (s  +  t)  in  the  neighborhood  of  the  origin 


